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Chapter  1 


Introduction 


The  precise  and  objective  measurement  of  the  severity  of  coronary  artery  obstruc¬ 
tions  is  important  in  the  treatment  of  patients  with  ischemic  heart  disease,  that  is, 
deficiency  of  blood  supply  to  the  heart. 

The  coronary  arteries  are  the  vessels  which  encircle  the  heart  and  supply  blood 
directly  to  the  heart  muscle.  Narrowing  of  the  coronary  arteries,  usually  as  a  result 
of  atherosclerosis,  prevents  adequate  supply  of  blood  to  the  heart,  causing  pain  and, 
eventually,  damage  to  the  heart  muscle.  Such  damage  of  the  heart  is  called  coronary 
disease.  The  number  of  deaths  from  coronary  disease  has  increased  strikingly  in  the 
last  70  years  in  association  with  a  marked  prolongation  of  the  average  life  span  [  1] . 
In  order  to  restore  normal  flow  of  blood  to  the  heart,  over  100,000  coronary  artery 
bypass  surgical  procedures  are  undertaken  each  year  [l]. 

An  angiogram  is  an  x-ray  picture  of  arteries  taken  after  a  contrast  agent  has 
been  injected  into  the  vessels.  Coronary  angiography  is  necessary  to  determine 
whether  a  bypass  operation  is  needed.  In  addition,  it  is  an  essential  part  of  any 
surgical  therapeutic  approach,  as  it  is  used  for  the  preparation  and  evaluation  of 
the  results  of  surgery.  An  example  of  a  coronary  angiogram  is  shown  in  figure  1.1. 
In  the  case  of  coronary  angiograms  the  injection  of  the  contrast  agent  is  done  via 
a  catheter  directly  to  the  arteries  to  achieve  high  concentrations.  At  present  the 
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best  storage  medium  for  the  coronary  angiograms  is  35mm  cinefilm,  because  of  its 
inherent  high-  resolution. 

Conventional  visual  evaluation  of  the  severity  of  coronary  obstructions  has  been 
characterized  by  large  inter-  and  intra-observer  variability  12, 36, 11, 24 i.  Thus,  the 
development  of  digital  processing  techniques  for  the  precise  and  objective  measure¬ 
ment  of  the  severity  of  coronary  obstructions  from  coronary  angiograms  is  important 
in  the  treatment  of  patients  with  ischemic  heart  disease. 

Coronary  artery  imaging  presents  special  problems,  because  of  their  location 
on  the  beating  heart  and  their  shape  and  size.  For  example,  intravenous  digital 
radiography  can  be  used  to  visualize  most  of  the  major  arteries  except  the  coro¬ 
naries.  When  the  injection  is  intravenous  the  SNR  of  the  angiographic  image  is 
low.  To  improve  the  SNR,  an  x-ray  image  obtained  before  injection  of  contrast 
agent  into  the  vessels  is  subtracted  from  an  image  obtained  after  injection.  The 
image  that  is  obtained  in  this  way  contains  the  vessels  in  uniform  background.  In 
the  coronary  artery  case  subtraction  results  in  poor  quality  images,  because  of  the 
rapid  movement  of  the  heart.  As  a  result  of  these  special  problems  many  techniques 
that  have  been  used  quite  successfully  for  stenosis  determination  of  other  arteries, 
like  the  femoral  and  carotid,  do  not  have  satisfactory  performance  when  applied  to 
coronary  arteries. 

The  most  commonly  used  measure  of  an  obstruction  is  percent  stenosis,  which 
is  determined  as  the  maximum  percent  arterial  narrowing  within  a  specified  length 
of  the  vessel 

stenosis  (%)  =  — — -  100  (1.1) 

where  d\  is  the  diameter  (or  cross-sectional  area)  of  the  normal  section  of  the  artery 
and  d2  is  the  diameter  (or  cross-sectional  area)  of  the  most  obstructed  section.  Of 
course,  the  resulting  stenosis  is  different  depending  on  whether  we  use  the  diameter 
or  the  cross-sectional  area.  According  to  equation  (1.1),  in  determining  the  severity 
of  an  obstruction,  it  is  the  relative  dimensions  of  the  vessel  that  matter,  not  the  ab- 


solute  ones.  However,  in  some  cases  it  is  important  to  have  absolute  measurements 
of  vessel  dimensions,  especially  when  it  is  not  possible  to  locate  with  certainty  a 
normal  section  of  the  vessel. 

The  hemodynamic  resistance  of  the  vessel  is  determined  by  the  cross-sectional 
area  of  the  vessel.  According  to  Poisenille’s  equation,  the  pressure  gradient  P 
necessary  to  ensure  a  constant  flow  F  in  a  cylindrical  vessel  of  radius  r  and  length 
L  is  given  bv 

P  =  ^  ^  (1-2) 

nr*  A- 

where  ^  is  the  coefficient  of  viscosity  and  A  is  the  cross-sectional  area  of  the  vessel  [4]. 
This  theoretical  result  has  also  been  corroborated  by  experimental  studies  [13,14]. 
Thus,  the  vessel  dimension  that  is  most  significant  for  this  problem  is  the  cross- 
sectional  area.  The  severity  of  heart  disease  is  determined  by  the  cross-sectional 
area  of  the  obstructed  segments  of  the  vessel. 

According  to  almost  unanimous  opinion,  the  stenoses  regarded  as  significant 
enough  to  require  an  operation  are  those  which  reduce  the  lumen  (that  is  the  cross- 
sectional  area)  by  7o°7  100°r  4  .  Thus,  it  is  important  to  have  accurate  estimates 
of  the  vessel  dimensions  at  its  most  narrow  points,  as  they  may  determine  whether 
an  operation  is  needed. 

The  cross-sectional  area  of  the  vessel  can  be  obtained  from  the  diameter,  if  the 
vessel  has  a  circular  cross-section.  This  is  approximately  true  for  healthy  vessels, 
but  at  a  coronary  obstruction  the  shape  of  the  cross-section  is  frequently  irregular. 
Figure  1.2  shows  a  few  characteristic  examples  of  such  cross-sections  taken  from 
Abrams  1  .  In  such  cases  it  is  obvious  that  the  diameter  depends  on  the  angle  of 
projection,  and  does  not  provide  a  good  estimate  of  the  area. 

The  methods  which  estimate  the  area  of  the  cross-sections  directly  from  the 
image  density  are  usually  called  densitometric.  These  methods  assume  that  the 
contrast  agent  is  mixed  homogeneously  with  blood,  and  also  that  there  is  a  linear 
relationship  between  the  density  of  the  radiographic  image  and  the  thickness  of  con- 


trast  agent.  This  relationship  can  be  secured  under  rigidly  controlled  experimental 
radiographic  “Conditions  (23.9i.  When  this  is  not  possible,  appropriate  corrections 
must  be  taken.  Spears  developed  a  rotating  wedge  technique  for  linearization  of 
cineangiographic  gray  scale  with  contrast  medium  thickness  30,31  .  A  different 
technique  for  transforming  x-ray  density  levels  to  irradiated  object  thickness  was 
implemented  by  Reiber  et  a I.  (17).  Thus,  densitometric  methods  require  strict 
densitometric  checks  and,  if  necessary,  calibration  of  the  imaging  system. 

The  maximum  percent  area  reduction,  as  defined  in  equation  (1.1),  is  a  reason¬ 
able  estimate  of  heart  disease.  Obviously,  a  more  accurate  indication  of  the  severity 
of  the  obstructions  would  be  the  exact  shape  of  the  cross-section  at  each  point  along 
the  artery.  As  we  will  see  in  chapter  2  the  optical  density  distribution  of  the  ra¬ 
diographic  image  is  the  line-integral  projection  of  the  x-ray  attenuation  coefficient 
along  each  ray  path  through  the  human  body.  According  to  the  projection-slice  the¬ 
orem,  if  we  know  the  projections  at  all  angles,  then  we  can  reconstruct  the  object 
exactly.  Since  the  attenuation  coefficient  of  the  contrast  agent  inside  the  arteries 
is  different  from  the  coefficient  of  the  surrounding  material,  we  can  thus  obtain 
the  vessel  cross-sections.  However,  there  are  practical  problems  in  obtaining  and 
combining  angiograms  of  more  than  one  view.  More  than  one  view  perpendicular 
to  the  long  axis  of  a  given  vessel  segment  often  cannot  be  obtained,  because  of  over¬ 
lapping  branches  31.33  .  Even  when  we  can  get  more  views,  proper  registration  of 
the  projections  in  the  various  views  is  difficult  to  achieve  because  of  motion  of  the 
heart.  Thus  it  is  important  to  extract  as  much  information  as  possible  about  the 
vessel  cross-sections  from  one  projection. 

It  is  clear  that  we  cannot  reconstruct  the  vessel  cross-sections  from  one  projec¬ 
tion.  Even  if  we  could  isolate  the  vessel  component  of  the  projection,  we  could  not 
reconstruct  the  vessel  cross-section,  unless  we  had  some  a  priori  information  about 
its  shape.  However,  some  vessel  parameters  can  easily  be  obtained  from  the  vessel 
component  of  the  projection.  Figure  1.3  shows  one  of  the  vessel  cross-sections  of 
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Figure  1.3:  Vessel  Cross-Section  and  Projection 

figure  1.2  and  its  projection,  assuming  a  perfect  imaging  system  and  uniform  con¬ 
trast  agent  concentration  within  the  vessel.  This  projection  corresponds  to  a  linear 
profile  perpendicular  to  the  two-dimensional  vessel  projection.  The  cross-sectional 
area  is  the  integral  along  this  profile.  The  diameter  of  the  cross-section  is  simply 
the  length  of  the  nonzero  projection  values.  It  is  clear  that  the  diameter  changes 
when  we  rotate  the  vessel,  while  the  cross-sectional  area  does  not. 

The  problem  of  course  is  that  we  do  not  have  the  vessel  component  of  the  pro¬ 
jection.  There  are  other  objects  in  the  background,  and  distortions  are  introduced 
by  the  imaging  system.  The  task  of  isolating  the  vessel  component  is  not  easy.  As  a 
result,  the  problems  of  estimating  the  diameter  and  cross-sectional  area  of  the  vessel 
become  nontrivial.  Figure  1.4(a)  shows  the  projection  of  figure  1.3,  which  assumes 
a  perfect  imaging  system.  Figure  1.4(b)  shows  the  same  projection  after  we  added  a 
low  order  polynomial  (background),  convolved  with  a  Gaussian  line  spread  function 
(blur),  and  added  white  noise.  The  resulting  profile  is  close  to  what  we  would  ob¬ 
tain  from  the  actual  imaging  system.  We  observe  that,  while  it  is  straight  forward 


(a)  Perfect  Imaging  System 


(b)  Actual  Imaging  System 


Figure  1.4:  Projection  of  Vessel  Cross-Section 

to  obtain  the  diameter  and  the  area  from  the  first  projection,  it  is  not  obvious  how 
we  would  obtain  the  same  parameters  from  the  second. 

1.1  Problem  Statement 

The  problem  we  will  consider  in  this  thesis  is  the  following.  Given  a  coronary 
angiogram,  that  is,  one  view  of  an  artery,  estimate  the  diameter  and  cross-sectional 
area  at  each  point  along  the  artery.  We  will  assume  that  we  have  a  single  isolated 
vessel.  We  will  not  attempt  to  handle  branch  points  and  overlapping  vessels.  Our 
procedure  for  computer  measurement  of  coronary  angiograms  assumes  that  the 
angiographer  selects  the  vessel  segment  to  be  analyzed. 

The  diameter  we  will  be  estimating  is  the  one  parallel  to  the  projection  plane.  As 
we  saw  above  this  diameter  is  sensitive  to  vessel  rotation.  The  cross-sectional  area 
on  the  other  hand  does  not  present  this  problem.  In  addition  to  this,  the  cross- 
sectional  area  determines  the  hemodynamic  resistance  of  the  vessel.  Hence  the 


cross-sectional  area  is  a  more  attractive  measure  of  vessel  obstruction.  However, 
cross-sectional -area  estimation  assumes  that  there  is  a  linear  relationship  between 
the  density  of  the  radiographic  image  and  the  thickness  of  contrast  agent.  As 
we  have  already  mentioned,  when  this  is  not  true,  appropriate  corrections  must 
be  taken.  The  advantage  of  the  diameter  estimates  is  that  they  are  not  sensitive 
to  deviations  from  this  linear  relationship.  The  experiments  with  vessel  phantom 
angiograms  in  chapters  4  and  5  illustrate  these  points. 

The  severity  of  heart  disease  is  determined  by  the  most  obstructed  segments  of 
the  vessel.  It  is  usually  the  greatest  degree  of  narrowing  observed  that  is  recorded 
by  the  angiographer  [4].  Thus  it  is  important  to  have  accurate  estimates  of  the 
vessel  dimensions  at  the  most  narrow  points.  An  algorithm  that  estimates  the 
vessel  dimensions  very  accurately  at  most  points  along  the  vessel  except  the  very 
narrow  ones  is  not  very  useful. 

1.2  Existing  Techniques 

Current  algorithms  analyze  the  angiograms  by  examining  the  densitometric  profiles 
of  a  series  of  lines  perpendicular  to  the  vessel  image.  Their  objective  in  most  cases 
is  to  find  the  diameter  or  the  area  of  the  vessel.  In  both  cases  the  boundary  points 
of  the  vessel  must  be  determined.  Most  algorithms  are  quite  heuristic  and  find 
the  boundaries  of  the  artery  by  first  smoothing  the  profile  and  then  finding  the 
inflection  points  (extrema  of  the  first  derivative),  or  the  knee  points  (maxima  of 
the  second  derivative)  [26,16,3,32,15],  Smoothing  is  done  by  averaging  or  fitting 
some  polynomial  or  spline  function  to  the  projection  profile.  Barth  et  a /.  [5j  use 
a  correlation  technique  to  find  the  positions  of  maximal  second  derivative.  Smith 
et  a/.  29  do  not  use  derivative  operators  to  determine  the  boundaries;  instead 

they  select  a  threshold  level  after  background  subtraction.  The  diameter  of  the 
vessel  is  obtained  as  the  distance  between  the  two  boundary  points.  The  area 


is  obtained  as  the  integral  of  the  projection  profile  between  the  boundary  points 
after  the  background  has  been  estimated  and  subtracted.  The  estimation  of  the 
background  is,  in  our  view,  a  nontrivial  task.  Some  authors  [16]  estimate  the 
background  by  computing  the  linear  regression  line  through  the  points  directly  left 
and  right  of  the  detected  boundaries.  Smith  et  a I.  29  estimate  the  background  by 
sliding  a  large  sphere  underneath  the  surface  of  the  image.  An  extensive  overview 
of  many  <  f  these  heuristic  methods  is  given  in  201. 

We  will  refer  to  the  methods  described  above  as  the  slope  methods,  because 
they  base  their  analysis  on  the  slope  of  the  profiles.  Their  main  advantage  is  their 
computational  simplicity.  They  are  heuristic,  however,  and  their  performance  is  not 
satisfactory.  As  we  will  see  in  chapter  4  this  is  especially  true  at  the  narrow  sections 
of  the  vessels,  which  are  the  most  important  ones.  We  will  also  show  that  they  are 
sensitive  to  scaling  of  the  image.  The  choice  of  the  inflection  points  is  arbitrary, 
except  for  the  fact  that  they  are  close  to  what  our  eyes  perceive  as  the  vessel  edge. 
In  fact,  the  inflection  points  are  typically  located  inside  the  actual  vessel  edges. 
The  knee  points  happen  to  be  closer  to  the  true  boundaries,  but  again  the  choice  is 
arbitrary.  As  one  would  expect,  the  knee  points  are  less  reliable  than  the  inflection 
points,  because  second  derivatives  are  involved. 

Another  problem  with  the  slope  methods  is  that  they  do  not  take  fully  into 
account  the  spatial  continuity  of  the  vessel.  One  way  to  use  the  vessel  continuity 
is  by  averaging  the  profiles  before  finding  the  edge  points.  An  alternative  is  to 
smooth  the  edge  points  after  detection.  Finally,  the  detected  points  of  one  profile 
can  be  used  to  determine  bounds  for  the  location  of  the  edge  points  of  the  next 
profile.  There  are  problems  associated  with  each  approach.  Averaging  of  the  profiles 
may  result  in  loss  of  resolution.  The  same  is  true  for  smoothing  of  the  detected 
boundaries.  The  use  of  the  bounds  can  lead  to  very  distorted  boundaries,  or  totally 
wrong  ones,  because,  if  the  algorithm  starts  tracking  the  wrong  edge,  it  is  impossible 
for  it  to  recover. 


Inaccuracies  in  the  estimation  of  the  vessel  boundaries  by  the  slope  methods 
affect  the  crOsS-sectional  area  estimates  as  well  as  the  diameter  estimates.  The 
area  estimates  are  affected  because  the  area  is  obtained  as  the  integral  between  the 
boundary  points  after  background  subtraction,  and  also  because  the  background 
estimate  depends  on  the  boundary  points. 

Shmueli  used  a  modeling  approach  to  detect  vessel  boundaries  in  subtraction 
angiograms  27,28  .  Subtraction  angiograms  are  obtained  when  an  x-ray  image 
obtained  before  injection  of  contrast  agent  into  the  vessels  is  subtracted  from  an 
image  obtained  after  injection.  The  purpose  of  the  subtraction  is  to  remove  the 
background  and  thus  increase  the  signal  to  noise  ratio  (SNR)  of  the  vessel  image. 
In  subtraction  angiography  the  injection  is  intravenous  and  so  the  SNR  of  the 
final  image  is  low.  Shmueli’s  algorithm  is  designed  to  detect  vessels  in  uniform 
background  and  was  applied  to  simulated  data  and  to  vessel  phantom  measurement 
data.  As  we  have  already  mentioned  subtraction  angiograms  of  the  coronary  arteries 
are  especially  difficult  to  obtain  and  are  not  currently  used  for  diagnostic  purposes, 
because  of  their  poor  quality  (low  SNR). 

Shmueli  assumes  a  simple  model  for  the  profiles  perpendicular  to  the  vessel.  It 
is  based  on  the  assumption  that  the  vessel  has  a  circular  cross-section  and  there  is 
no  background.  Thus  the  profiles  are  characterized  by  the  radius  and  the  location 
of  the  center.  Shmueli  uses  a  stochastic  approach  to  model  vessel  continuity. 

Some  authors  have  developed  algorithms  for  the  reconstruction  of  the  vessel 
cross-sections  from  two  orthogonal  projections  [16;.  They  use  some  heuristic  tech¬ 
nique  to  find  the  vessel  boundaries  and  to  remove  the  background  from  each  pro¬ 
jection,  and  then  they  reconstruct  the  cross-sections  in  the  following  way.  They 
assume  that  the  cross-section  can  be  represented  by  a  matrix  of  0’s  and  l’s,  and 
that  the  projections  are  proportional  to  the  row  and  column  sums  of  the  matrix. 
Since  it  is  well  known  that  the  matrix  can  not  be  reconstructed  uniquely  from  the 
two  sums,  their  reconstruction  algorithm  minimizes  a  cost  criterion  that  reflects  the 


fact  that  adjacent  slices  must  be  similar. 

The  probTem  with  this  approach  is  that,  in  addition  to  the  fact  that  the  cross- 
sections  cannot  be  reconstructed  from  the  two  projections,  the  algorithm  assumes 
that  the  vessel  projections  can  be  isolated  from  the  background.  We  have  already 
discussed  the  problems  with  the  estimation  of  boundaries  and  background. 


1.3  Scope  of  the  Thesis 

In  this  thesis  we  present  a  modeling  approach  to  estimate  the  dimensions  of  the 
coronary  arteries  from  their  x-ray  projections.  The  idea  is  the  following.  We  fit 
a  model  to  the  x-ray  image  and  obtain  the  desired  vessel  dimensions  from  the 
estimated  model  parameters,  instead  of  relying  on  direct  boundary  estimates  based 
on  local  extrema  of  first  or  second  derivatives. 

We  develop  a  two-dimensional  model  of  a  rectangular  section  of  the  angiogram 
that  contains  a  single  vessel.  Our  model  accounts  for  the  structure  of  the  back¬ 
ground,  as  well  as  distortions  introduced  by  the  imaging  system.  The  desired  vessel 
dimensions,  that  is,  the  diameter  and  cross-sectional  area  at  each  point  along  the 
vessel,  can  be  obtained  easily  from  the  model  parameters. 

The  vessel  component  of  the  projection  is  derived  from  a  simple  three  dimen¬ 
sional  model  of  the  vessel.  The  vessel  is  assumed  to  be  a  tube  with  elliptical  cross- 
sections.  The  background  is  modeled  by  a  low  order  two-dimensional  polynomial. 
The  model  takes  into  account  the  continuity  of  both  the  vessel  and  the  background. 
In  addition  to  the  two-dimensional  model  of  a  section  of  the  angiogram,  we  derive 
a  one-dimensional  model  of  the  densitometric  profiles  of  lines  crossing  the  vessel. 
This  model  can  also  be  used  to  obtain  the  diameter  and  cross-sectional  area  at  each 
point  along  the  vessel,  even  though  it  does  not  use  the  continuity  of  the  vessel  and 
the  background. 

The  main  difference  with  Shmueli’s  model  is  the  addition  of  the  background. 


T 


Shmueli’s  model  assumes  a  flat  background  and  therefore  cannot  even  be  applied  to 
coronary  angiograms.  Also  Shmueli's  model  of  vessel  continuity  is  stochastic.  The 
model  presented  in  this  thesis  is  deterministic. 

The  parameters  of  the  model  are  deterministic  and  we  obtain  the  maximum 
likelihood  estimates,  which  is  equivalent  to  finding  the  model  that  best  fits  the  data 
in  the  least  squares  sense.  The  estimation  of  the  model  parameters  is  a  non-linear 
problem.  We  developed  an  iterative  algorithm  to  estimate  the  parameters.  The 
computational  requirements  for  the  one-dimensional  algorithm  are  reasonable,  even 
though  they  are  considerably  greater  than  the  slope  methods.  The  two-dimensional 
algorithm  requires  much  more  computation  than  the  one-dimensional  approach. 

We  test  the  algorithms  on  three  types  of  data.  The  first  type  is  computer  gener¬ 
ated  data.  It  provides  the  opportunity  to  test  the  algorithms  under  the  assumption 
that  the  model  is  valid  and  to  compare  them  with  the  current  methods.  Second, 
the  algorithms  are  tested  on  x-rays  of  contrast-medium-filled  cylindrical  phantoms. 
Since  the  true  vessel  dimensions  are  known,  this  set  of  data  gives  us  the  opportunity 
to  test  the  model  of  the  imaging  system.  Finally,  the  algorithms  are  tested  on  real 
coronary  angiograms.  In  this  case  the  true  vessel  dimensions  are  not  known.  Thus 
this  set  of  experiments  can  only  lead  to  qualitative  conclusions.  In  the  first  part  of 
our  experiments  we  compare  the  performance  of  this  one-dimensional  model  to  the 
slope  approaches.  We  then  compare  the  performance  of  the  one-dimensional  model 
with  that  of  the  two-dimensional  model. 

Our  results  indicate  that  both  the  one-dimensional  and  the  two-dimensional 
algorithm  have  better  performance  than  the  conventional  slope  approaches.  This  is 
because  the  modeling  approach  uses  all  the  information  in  the  projection  image  to 
estimate  the  model  parameters.  The  slope  approaches  use  only  the  slope  information 
at  the  sides  of  the  vessel  to  determine  the  vessel  boundaries  and  get  an  estimate  of 
the  background.  We  also  show  that  the  two-dimensional  algorithm  is  more  robust 
to  non-white  noise  in  the  background  than  the  one-dimensional  algorithm. 


1.4  Summary  and  Notation 

The  thesis  is  organized  as  follows.  The  parametric  model  approach  for  the  detection 
of  coronary  artery  dimensions  is  presented  in  chapter  2.  The  computation  of  the 
model  parameters  is  discussed  in  chapter  3.  In  chapter  4  we  compare  the  modeling 
approach  to  the  slope  methods.  The  two  versions  of  the  modeling  approach  are 
compared  in  chapter  5.  The  conclusions  of  the  thesis  are  summarized  in  chapter  6. 

Throughout  the  thesis  we  will  use  parentheses  to  denote  a  continuous  function 
/(.)  and  brackets  to  denote  a  discrete  function  /j.j.  We  will  use  bold  face  lower  case 
letters  to  denote  a  vector  x. 


Chapter  2 


Coronary  Angiogram  Model 


The  radiographic  imaging  chain  includes  the  following  stages,  shown  in  figure  2.1. 
The  object  is  illuminated  by  the  x-ray  source  into  the  image  intensifier.  A  cine¬ 
camera  records  the  output  of  the  image  intensifier  on  35mm  film.  Each  frame  is 
optically  magnified  and  digitized  by  a  vidicon  camera/digitizer. 

First  we  will  develop  a  model  for  the  coronary  angiogram  under  the  assumption 
of  a  perfect  imaging  system.  That  is,  the  x-ray  tube  is  an  ideal  point  source  emitting 
a  monoenergetic  x-ray  beam,  the  image  intensifier  introduces  no  distortions,  the 
film  characteristic  curve  is  linear,  there  is  no  film  grain  noise,  and  the  output  of 
the  system  is  a  continuous  image.  Later  we  will  extend  our  model  to  take  into 
consideration  the  more  significant  distortions  of  the  imaging  system. 


Figure  2.1:  Imaging  System 


2.1  Perfect  Imaging  System 
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The  analysis  of  the  x-ray  images  is  based  on  the  Lambert-Beer  law  [23].  According 
to  this  law,  a  monoenergetic  x-ray  beam  with  input  intensity  h  photons  passing 
through  a  length  2  of  material  gives  an  output  intensity  of 

/  =  he (2.1) 

The  attenuation  coefficient  /x  depends  on  the  energy  of  the  x-ray  beam  and  the 
type  of  material  it  passes  through.  In  the  case  of  the  human  body,  the  attenuation 
coefficient  /x(x,y,2)  is  a  function  of  position,  and  the  simple  fiz  product  becomes  a 
line  integral.  Thus  the  output  intensity  is  also  a  function  of  position 

/(i,y)  =  /o(x,y)e-/^^d*  (2.2) 

The  projection  P(x,y)  that  we  see  in  an  x-ray  image  is  actually  the  following  quan¬ 
tity  25] 

P(x,  y)  =  -  In  y  -'  =  [  n(x,  y,  2)  dz  (2.3) 

/o(x,y)  J 

Thus  the  x-ray  image  is  the  line-integral  projection  of  the  attenuation  coefficients. 
The  x-ray  beam  passing  through  the  human  body  encounters  three  types  of  material: 
tissue,  bones,  and  contrast  agent  (iodine)  which  has  been  injected  into  the  blood 
vessels.  Thus 


P{x,y)  =  j^i{x,y,z)dz  -t-  J^nT(x,y,  z)  dz  +  j^y.B(x,y,z)  dz  (2.4) 


where  p/(x,y,2)  is  the  attenuation  coefficient  of  iodine,  and  the  corresponding  in¬ 
tegral  is  over  the  iodine  sections  along  the  line  parallel  to  the  z-axis  at  (x,y).  The 
other  two  integrals  are  similarly  defined.  If  we  assume  that  iodine  is  uniformly 
distributed  inside  the  vessel,  and  thus  the  attenuation  coefficient  is  constant  in  the 
iodine  sections,  then  we  have 


P(x,y)  =  M/0(x,y)  +  J^nb(x,y,z)dz 


(2.5) 
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where  //(x,  y)  is  the  thickness  of  the  blood  vessel  (iodine)  along  the  line  parallel  to 
the  z-axis  at"(j f,  y)  of  the  body,  and  the  integral  is  the  combination  of  the  tissue  and 
bone  terms  of  equation  (2.4).  Thus  the  projection  consists  of  two  terms,  the  vessel 
term  and  the  background  term,  and  equation  (2.5)  can  be  written  as 

P(i,y)  =  v(i,y)  -  6(x,y)  (2.6) 

The  attenuation  coefficient  of  the  contrast  agent  is  substantially  higher  than  the 
other  two,  and  that  is  why  we  can  see  the  vessels. 

Assuming  this  perfect  image  processing  system  we  will  develop  a  model  for  the 
x-ray  image.  In  particular,  we  will  develop  a  model  for  the  vessel  component  and 
a  model  for  the  background  component.  The  background  model  will  be  for  a  small 
rectangular  section  of  the  image  around  the  vessel. 

2.2  Two-Dimensional  Vessel  Model 

In  this  section  we  will  develop  a  two-dimensional  model  for  the  vessel  component 
v(i,y)  of  equation  (2.6).  We  will  do  this  by  considering  a  model  for  the  three- 
dimensional  vessel  and  from  this  deriving  the  model  of  its  two-dimensional  projec¬ 
tion  on  the  x-ray  film. 

We  consider  the  following  model  for  the  three-dimensional  vessel.  We  assume 
that  the  vessel  is  continuous  and  has  elliptical  cross-sections.  Then  it  can  be  de¬ 
scribed  by  a  continuous  multi-dimensional  function,  each  point  of  which  is  a  vector 
containing  the  coordinates  of  the  axis  and  the  parameters  characterizing  the  el¬ 
liptical  cross-section  perpendicular  to  the  axis  at  the  point.  This  kind  of  object 
is  usually  referred  to  as  a  generalized  cone  or  cylinder  {2  .  Figure  2.2  shows  an 
example  of  such  a  cylinder  and  its  line-integral  projection. 

This  model  is  clearly  too  simple  to  match  the  actual  vessel  exactly.  The  model  is 
pretty  close  to  the  healthy  vessels,  but  it  deviates  significantly  from  the  obstructed 
vessels.  As  we  pointed  out  in  the  introduction,  the  shape  of  a  cross-section  at  an 


Figure  2.2:  Generalized  Cylinder  and  Its  Projection 


obstruction  cai.  be  irregular  in  shape.  The  model  is  simple  so  that  we  can  estimate 
its  parameters  from  only  one  projection.  The  hope  is  that  even  when  the  model 
is  crude,  we  can  still  get  good  estimates  of  the  true  vessel  parameters,  such  as  the 
diameter  and  cross-sectional  area. 

A  slice  of  a  cylinder  with  an  elliptical  cross-section  is  always  an  ellipse.  Thus, 
if  the  vessel  is  continuous  and  its  parameters  are  slowly  varying,  then  the  slices  of 
the  vessel  are  ellipses.  Now.  look  at  a  plane  perpendicular  to  the  x-ray  image  which 
slices  the  vessel,  and  consider  the  line-integral  projection  of  the  vessel  along  that 
plane.  The  elliptica1  vessel  cross-section  and  its  projection  are  shown  in  figure  2.3 
in  solid  lines. 

An  ellipse  is  described  by  the  following  parameters:  the  major  and  minor  axes, 
the  angle  of  rotation,  and  the  coordinates  of  its  center.  If  the  axes  have  lengths  2a 
and  26,  the  angle  of  rotation  is  0 .  and  the  center  of  the  ellipse  is  located  at  (c(,cs), 
then  the  equation  of  the  ellipse  is 

(.s  -  c,)  cos$  -  (t  -  c()  sind  2  —  (s  —  c.)  stnd  ~  (t  -  c,)  cosd  2 
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and  its  projection  is  given  by 


p(()  =  hrvrI  «t-c,\<r  (28) 

0,  otherwise 

where  r2  =  b2cos26  —  a2sin26. 

It  is  easy  to  see  that  this  ellipse  has  the  same  projection  as  an  ellipse  with  axes 
a  =  ab/r  and  b  =  r  and  6  =  0.  The  new  ellipse  is  shown  in  figure  2.3  in  dotted  line. 
This  shows  that,  as  we  mentioned  earlier,  it  is  impossible  to  recover  the  exact  shape 
of  the  cross-section  from  one  projection.  However,  the  area  of  the  two  ellipses  is  the 
same  and  is  given  by 

area  =  nab  =  nab  (2.9) 


which  is  also  equal  to  the  integral  of  the  projection. 

In  fact,  if  we  express  the  projection  of  an  ellipse  in  terms  of  the  parameters  r, 
a.  and  c  as  follows, 


v(0  = 


2 at\Jr2  -  (t  -  c)2,  if  t  -  cj  <  r 
0,  otherwise 


(2.10) 


then  these  parameters  can  be  recovered  from  the  projection.  The  area  of  the  cross- 
section  is  given  by 

area  =  7rar2  (2.11) 

Thus  two  quantities  that  we  can  get  from  one  projection  are  the  area  of  the  cross- 
section,  and  its  diameter  2r,  which  is  equal  to  the  diameter  of  the  ellipse  which  has 
the  given  projection  and  an  axis  parallel  to  the  projection  plane. 

As  an  illustration  of  the  potential  errors  in  the  obstructed  sections  of  the  vessels 
we  considered  several  cross-sections  obtained  from  post-mortem  analyses.  Figure  2.4 
shows  the  line-integral  projections  of  the  cross-sections  of  figure  1.2.  Also  shown  in 
this  figure  are  the  elliptical  models  fitted  to  these  projections.  The  model  parameters 
are  chosen  so  that  the  sum  of  squares  of  the  difference  between  the  actual  projection 
and  the  model  is  minimized.  Table  2.1  displays  the  percent  error  between  the 


estimated  and  actual  parameters.  We  see  that,  even  though  the  cross-sections  are 
far  from  being,  elliptical,  we  get  good  estimates  of  the  cross-sectional  area.  The 
diameter  estimates  are  not  as  good,  but  the  exact  values  are  not  very  useful,  since 
they  depend  on  the  angle  of  projection. 

So  far  we  have  derived  a  model  for  a  cross-section  of  the  vessel  projection.  We 
must  now  extend  this  to  a  two-dimensional  model  for  the  projection  of  the  vessel 
on  the  x-ray  image.  Consider  the  profiles  of  the  projection  that  are  perpendicular 
to  the  axis  of  the  projected  vessel.  According  to  our  discussion,  each  profile  is 
described  by  an  equation  of  the  form  of  (2.10).  Let  ci {q)  and  c2(q)  be  the  functions 
that  specify  the  x-  and  y-coordinates  of  the  axis  of  the  projected  vessel,  and  let  the 
functions  r(q)  and  a(q)  describe  the  corresponding  perpendicular  profiles.  These 
functions  describe  the  vessel  component  v(x,y)  of  the  perfect  projection  model  of 
equation  (2.6).  The  continuity  of  the  vessel  implies  that  the  functions  r(<j),  a(q), 
cy(q)  and  c2(<f)  are  continuous.  We  will  assume  that  they  are  cubic  splines. 

Splines  are  piecewise  polynomial  functions  of  degree  n  that  are  connected  to¬ 
gether  (at  points  called  nodes)  so  as  to  have  n  -  1  continuous  derivatives  [8,10]. 
They  represent  nice  smooth  curves  which  are  convenient  to  use  in  many  practical 
applications.  One  important  property  of  splines  is  that  they  are  good  approxima¬ 
tions  of  functions  whose  behavior  in  one  region  is  totally  unrelated  to  their  behavior 
in  another  region.  The  behavior  of  polynomials,  as  well  as  most  other  mathemati¬ 
cal  functions,  is  just  the  opposite.  Their  behavior  in  any  small  region  affects  their 
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behavior  everywhere  [21]. 

The  continuity  and  flexibility  of  the  cubic  splines  are  the  main  reasons  that  make 
them  a  good  choice  for  modeling  the  spatial  continuity  of  the  vessel.  In  addition, 
their  flexibility  can  be  adjusted  by  changing  the  spacing  of  the  nodes.  The  choice  of 
the  spline  parameters  does  not  depend  on  the  length  of  the  vessel.  The  importance 
of  these  properties  will  become  more  evident  in  the  next  chapter,  where  we  discuss 
the  estimation  of  the  model  parameters. 

We  showed  that  from  the  projection  parameters  we  can  obtain  the  area  of  the  ver¬ 
tical  cross-sections  of  the  three-dimensional  vessel,  but  not  the  exact  cross-sections. 
It  should  also  be  obvious  that  we  cannot  recover  the  exact  three-dimensional  loca¬ 
tion  of  the  vessel  either. 

Here  we  should  note  that  the  parameter  a  is  also  proportional  to  the  attenuation 
coefficient  p/  of  the  contrast  agent,  i.e.  to  the  iodine  concentration.  Assuming  that 
the  iodine  concentration  is  constant  along  the  vessel  section  of  interest,  estimates 
of  the  vessel  model  parameters  will  give  us  relative  estimates  of  the  true  vessel 
parameters. 


2.3  Two-Dimensional  Model  of  Background 

In  this  section  we  will  introduce  a  two-dimensional  model  for  the  background  com¬ 
ponent  b{x,y)  of  equation  (2.6).  We  will  model  the  background  in  a  rectangular 
section  of  the  image  around  the  vessel. 

We  assume  that  the  objects  in  the  background,  that  is  tissue  and  bones,  are 
much  bigger  than  the  blood  vessels.  Therefore  the  density  of  their  two-dimensional 
projection  is  slowly  varying,  and  we  can  model  it  by  a  low  order  two-dimensional 
polynomial  p(x,y).  A  two-dimensional  polynomial  is  a  function  of  the  form 

(2.12) 

i=0 ;=0 
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If  we  assume  that  a,  g_,  *  0  for  some  0  <  i  <  Q,  then  the  degree  of  the  polynomial 
is  Q  and  th$  prder  is  Q  1.  Given  that  the  image  occupies  a  rectangular  area 
a.b  ■  c,d  of  the  plane,  we  ran  express  this  function  in  terms  of  an  orthonormal 
basis  c\,(x.y) 

Q  Q-t 

P(x<y)  =  EH  dt30,j(x,y)  (2.13) 

i=o  j=0 

where  ipt}(x.y)  =  0*(x)<pv](y)  and  <t>f{x),  <t> *(y)  are  the  orthonormal  bases  for  the 
one-dimensional  polynomials  in  the  intervals  ja,6i  and  c,d]  respectively  8j. 

Polynomials,  especially  of  low  order,  are  very  smooth  functions.  They  are  unable 
to  take  on  sharp  bends  followed  by  relatively  flat  behavior.  Thus  their  properties 
are  very  different  from  those  of  the  vessel  projection,  which  varies  drastically  in 
the  direction  perpendicular  to  its  axis.  This  makes  the  models  of  the  background 
and  vessel  relatively  independent,  and  thus  makes  simultaneous  estimation  of  their 
parameters  from  one  projection  possible. 

The  model  for  the  perfect  projection  of  equation  (2.6)  now  becomes 

P(x,y)  =  v(x,y)  -  p(x,y)  (2.14) 

An  example  of  the  perfect  model  and  its  components  is  shown  in  figure  2.5. 
Figure  2.5(a)  depicts  the  vessel  model.  A  fifth  degree  polynomial  modeling  the 
background  is  shown  in  figure  2.5(b).  The  sum  of  the  two  components  is  shown  in 
figure  2.5(c). 

We  must  now  extend  our  model  to  take  into  consideration  the  distortions  of  the 
imaging  system. 

2.4  Model  of  Imaging  Process 

In  the  development  of  our  model  we  assume  that  the  x-ray  tube  is  an  ideal  point 
source  and  that  the  radiation  is  monoenergetic.  We  model  the  blurring  that  is 
introduced  at  various  stages  in  the  imaging  process  as  convolution  of  the  perfect 
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projection  image  with  a  two-dimensional  point  spread  function  y(x,y).  In  this  thesis 
we  will  assurae-that  the  point  spread  function  is  Gaussian  with  covariance  matrix  S , 
but  any  other  smooth  bell  shaped  function  can  be  used.  In  fact,  the  x-ray  pattern 
in  the  plane  of  the  film  (intensity  domain)  is  the  convolution  of  the  x-ray  pattern 
that  enters  the  image  intensifier  with  the  point  spread  function  of  the  film-screen 
combination  122:.  Tests  with  vessel  phantom  x-rays  show  that  it  is  reasonable  to 
model  this,  as  well  as  other  effects,  as  convolution  in  the  density  domain.  We  model 
degradations  such  as  quantum  mottle  noise,  digitization  noise,  and  film  grain  noise, 
as  white  noise  tu(x,y)  added  to  the  projection  image. 

According  to  this  model  the  optical  density  distribution  of  the  radiographic 
image  is  the  following 


f(x,y)  =  P(x,y)  *  g(x,y)  ~  w(x,y ) 


(2.15) 


where  P(x.y)  is  the  perfect  projection  of  equation  (2.14),  y(x,y)  is  the  Gaussian 
blurring  function,  and  w(x,y)  is  the  additive  noise.  Substituting  for  P(x.y)  we 
obtain  the  density  at  each  point  of  the  image  section 


/(x,  y)  =  [t>(x,  y)  4-  p(x,  y) i  *  y(x,  y)  +  u/(x,  y) 


(2.16) 


where 

/(x,y)  :  observed  densitometric  profile 
i '(x,y)  :  vessel  model 

p(x,y)  :  polynomial  model  of  background  tissue  and  bones 
y(x,y)  :  Gaussian  point  spread  function 
u'(x. y):  white  noise 

The  convolution  of  an  n-th  order  polynomial  with  a  smooth  function  like  a  Gaussian 
is  also  an  n-th  order  polynomial,  i.e. 


p(x,y)  *  g(x,y)  =  p'{x,y) 


(2.17) 


,y.y.V7i 


where  p'(x,y)  is  a  polynomial  of  the  same  order  as  p(x,y).  Substituting  equa¬ 
tion  (2.17)  into-(2.16)  we  have 

/( x,y)  =  v(x,y)  *  g(x,  y)  +  p'(x,  y)  +  u>(x,y)  (2.18) 

This  is  the  two-dimensional  model  of  the  coronary  angiogram  that  we  will  use 
in  this  thesis.  An  example  of  the  model  and  its  components  is  shown  in  figure  2.6. 
Figure  2.6(a)  depicts  the  vessel  model  of  figure  2.5(a)  convolved  with  a  Gaussian 
point  spread  function.  The  sum  of  this  component  and  the  background  of  fig¬ 
ure  2.5(b)  is  shown  in  figure  2.6(b).  When  white  noise  is  added  to  the  above  we 
obtain  figure  2.6(c). 

The  radiographic  imaging  process  is  much  more  complicated,  however.  There 
are  many  distortions  introduced  because  of  the  x-ray  tube  characteristics  (it  is 
not  an  ideal  point  source),  the  image  intensifier  distortions  (pincushion),  the  film 
characteristics,  and  the  fact  that  the  radiation  is  not  monoenergetic.  A  section  of 
an  actual  coronary  angiogram  is  shown  in  figure  2.7.  The  model  parameters  for 
figure  2.6  were  obtained  from  this  angiogram. 

Our  model  is  only  a  crude  model  and  cannot  account  for  all  the  distortions 
introduced  during  the  imaging  process.  However,  the  more  important  distortions 
are  either  accounted  for  by  the  model  (blurring,  noise),  or  can  be  corrected  inde¬ 
pendently.  Nonlii.oar  geometric  distortions,  such  as  the  pincushion  distortion,  are 
important  only  if  absolute  measurements  are  to  be  derived  from  the  angiograms. 
Kooijmar  et  a I.  15  present  a  technique  for  correcting  this  type  of  distortions.  Our 
model  assumes  a  linear  relationship  between  the  density  of  the  radiographic  image 
and  the  thickness  of  the  contrast  medium.  When  this  is  not  true,  linearity  can  be 
restored  by  the  techniques  we  mentioned  in  the  introduction  [30,31,171.  A  compar¬ 
ison  of  the  model  of  figure  2.6(c)  with  the  actual  angiogram  of  figure  2.7  indicates 
that  our  model  is  reasonable  despite  the  simplifications. 

As  we  discussed  in  the  introduction,  the  x-rays  of  contrast-medium-filled  cylin¬ 
drical  phantoms  give  us  an  opportunity  to  test  our  model.  Since  the  vessel  phantoms 
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Figure  2.7:  Real  Angiogram 


are  straight  with  circular  cross-sections,  and  the  background  is  flat,  they  both  fit  the 
model,  and  thus  we  can  use  their  x-rays  to  check  the  model  of  the  imaging  process. 
A  section  of  a  vessel  phantom  x-ray  is  shown  in  figure  2.8(a).  The  model  for  this 
phantom  is  shown  in  figure  2.8(b).  The  vessel  component  is  shown  in  figure  2.8(c), 
and  the  difference  between  the  original  and  the  vessel  component  is  shown  in  fig¬ 
ure  2.8(d).  Examining  this  figure  we  find  no  traces  of  the  vessel.  This  indicates 
that  the  model  of  the  imaging  process  is  adequate  for  this  case. 

Further  justification  for  our  model  will  be  provided  in  the  experiments  of  the 
following  chapters. 


(a)  Original 


(b)  Model 


"essel  Component 


(d)  Model  Minus  Vessel  Component 


Figure  2.8:  Comparison  of  Vessel  Phantom  and  Model 


2.5  One-Dimensional  Model 


r 


As  we  saw  in  the  introduction,  all  of  the  current  algorithms  estimate  the  vessel 
parameters  by  analyzing  the  densitometric  profiles  of  the  lines  perpendicular  to 
the  vessel.  A  model  of  these  one-dimensional  profiles  can  be  developed  based  on 
assumptions  similar  to  those  of  the  previous  sections.  Modeling  the  perpendicular 
profiles  would  be  an  obvious  first  step  towards  a  modeling  approach  to  our  problem. 
It  is  therefore  of  interest  to  compare  it  to  the  two-dimensional  model.  It  also  offers 
a  computational  advantage  over  the  two-dimensional  approach,  as  we  will  see  in  the 
next  chapter.  The  advantage  of  the  two-dimensional  approach  is  that  it  uses  the 
spatial  continuity  of  both  the  vessel  and  the  background. 

The  model  for  the  one-dimensional  profiles  can  be  derived  similarly  to  the  two- 
dimensional  model.  Consider  the  profile  of  a  line  crossing  the  vessel  projection  at 
some  angle,  not  necessarily  a  right  angle.  The  one-dimensional  analog  of  equa¬ 
tion  (2.16)  is  the  following 


/( 0  =  MO  -  p(0i  *  </(0  MO 


(2.19) 


where 


f(t)  :  observed  densitometric  profile 
v(t)  :  vessel  model 

p(0  :  polynomial  model  of  background  tissue  and  bones 
<7(f)  :  Gaussian  point  spread  function 
w(t):  white  noise 

The  model  for  a  cross-section  of  the  vessel  projection  v{t)  was  derived  in  section  2.2, 

and  is  given  by  equation  (2.10) 

.  ,  (  2 Q\/r2  -  (<  -  c)2,  if  t  -  c|  <  r 

f(0  =  v 

1 0,  otherwise 

The  background  is  modeled  as  a  low-order  polynomial  p(t)  of  the  form 

P(  0  =  XM'MO 


(2.20) 


where  <£,(<)  is  the  orthonormal  basis  for  the  one-dimensional  polynomials  in  an 
interval  a,6h-  -It  is  easy  to  see  that  a  slice  of  a  two-dimensional  polynomial  is  a 
one-dimensional  polynomial  of  the  same  degree.  Thus  the  vessel  and  background 
components  are  consistent  with  the  two-dimensional  model.  Blurring  is  modeled 
as  convolution  of  the  perfect  projection  v (t)  ~  p(t )  with  a  Gaussian  point  spread 
function  g(t ),  with  standard  deviation  s.  However,  it  is  not  true  that  a  slice  of  the 
convolution  of  two  functions  is  the  convolution  of  their  slices.  3!urring  in  all  direc¬ 
tions  cannot  be  equivalent  to  blurring  along  one  direction.  Thus  equation  (2.19)  is 
an  analog  of  equation  (2.16).  but  cannot  be  obtained  as  a  slice  of  equation  (2.16). 
The  final  component  of  the  model  is  additive  white  noise  w(t ),  with  standard  devi¬ 
ation  a. 

Similarly  to  the  two-dimensional  case  the  convolution  of  an  n-th  order  poly¬ 
nomial  with  a  smooth  function  like  a  Gaussian  is  also  an  n-th  order  polynomial, 
i.e. 

p(O*0(O=p'(O  (2-21) 

where  p'(t)  is  a  polynomial  of  the  same  order  as  p(t).  Substituting  equation  (2.21) 
into  (2.19)  we  have 

/(<)  =  v(t)  *  g(t)  -t  p'(t)  -  w(t)  (2.22) 

This  is  the  one-dimensional  model  of  the  densitometric  profile  along  a  line  crossing 
the  coronary  arter.,.  We  should  note  here  that  the  profile  does  not  have  to  be 
perpendicular  to  the  vessel  for  the  model  to  be  valid.  The  reason  for  which  we  are 
interested  in  the  perpendicular  profiles  is  that  it  is  their  parameters  that  are  useful 
for  our  problem. 

The  model  and  its  components  are  shown  in  figure  2.9.  Figure  2.9(a)  depicts 
the  vessel  model.  This  model  convolved  with  a  Gaussian  point  spread  function  is 
shown  in  figure  2.9(b),  while  a  fifth  degree  polynomial  modeling  the  background 
is  in  figure  2.9(c).  The  sum  of  the  components  of  figures  2.9(b)  and  2.9(c)  is  in 
figure  2.9(d).  When  white  noise  is  added  to  the  above  we  obtain  figure  2.9(e). 


An  example  of  an  actual  vessel  profile,  taken  from  a  real  angiogram,  is  shown  in 
figure  2.9(f).-  The  parameters  for  the  model  were  obtained  from  this  profile.  A 
comparison  of  the  model  of  figure  2.9(e)  with  the  actual  angiogram  of  figure  2.9(f) 
indicates  that  our  model  is  reasonable  despite  the  simplifications. 

As  we  did  in  the  two-dimensional  case,  we  can  test  the  model  of  the  imaging 
process  on  x-rays  of  contrast-medium-filled  cylindrical  phantoms.  A  profile  of  a  ves¬ 
sel  phantom  x-ray  is  shown  in  figure  2.10(a).  The  model  for  this  phantom  is  shown 
in  solid  line  in  figure  2.10(b);  the  original  is  shown  in  dotted  line  for  comparison. 
We  see  that  there  is  a  close  match  between  the  two.  The  vessel  component  is  shown 
in  figure  2.10(c),  and  figure  2.10(d)  shows  the  original  minus  the  vessel  component. 
No  trace  of  the  vessel  is  evident  in  this  subtracted  profile.  This  indicates  that  the 
model  of  the  imaging  process  is  adequate  for  this  case. 

The  two-dimensional  vessel  projections  can  be  described  by  the  parameter  func¬ 
tions  r(q)  and  a(q),  as  well  as  the  corresponding  axis  coordinate  functions  ci(g)  and 

The  difference  between  the  one-dimensional  and  the  two-dimensional  model  is 
that  the  latter  takes  into  consideration  the  continuity  of  both  the  vessel  and  the 
background.  Also  the  blurring  is  in  all  directions,  rather  than  along  a  linear  profile. 
In  the  next  chapter  we  will  consider  a  simple  heuristic  (and  suboptimal)  way  to 
incorporate  the  vessel  continuity  into  the  one-dimensional  model.  But  there  does 
not  appear  to  be  a  reasonable  way  to  use  the  continuity  of  the  background. 


Chapter  3 


Estimation  of  Model 


Parameters 


In  the  previous  chapter  we  introduced  parametric  models  for  a  rectangular  section 
of  the  coronary  angiogram  and  a  linear  profile  of  the  angiogram  crossing  the  artery. 
In  this  chapter  we  will  look  at  the  problem  of  estimating  the  parameters  of  these 
models  from  the  angiograms.  The  unknown  parameters  of  the  model  are  assumed 
to  be  nonrandom.  In  this  case,  the  problem  can  be  formulated  as  a  maximum 
likelihood  (ML)  estimation  problem. 

According  to  our  model  the  observations  (coronary  angiogram)  are  described  by 
equations  (2.18)  and  (2.22),  which  we  will  combine  into 


/(x)  =  /o(x;  a)  +  w{x) 


(3.1) 


where  /0(x;a)  =  v(x)  *<?(x)  +p'(x).  The  vector  x  stands  for  the  coordinate  pair  x, y 
in  the  two-dimensional  model,  or  the  coordinate  t  in  the  one-dimensional  model, 
and  a  is  the  vector  of  all  the  model  parameters.  The  parametric  function  /o  is 
deterministic,  and  the  noise  w  is  stationary,  white,  Gaussian,  with  variance  a2. 

In  the  Gaussian  case  the  ML  estimate  a  is  the  value  of  a  which  maximizes  the 
likelihood  function  1 34 1 


A[/(x); aj  =  j {/(x)  -  /0(x;a)}2dx 


(3.2) 
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Differentiating  with  respect  to  a,  we  find  that  the  ML  estimate  satisfies  the  equation 

I {/(x)  -  /0(x;a)}d— j*—  ^  _  Q  (3.3) 

The  covariance  matrix  of  the  ML  estimate  a  satisfies  the  Cramer-Rao  inequality 

E{(a-  a)(a  -  a)T}  >  J-1(a)  (3.4) 

where 


The  ML  equation  (3.3)  is  nonlinear  and  cannot  be  solved  analytically.  Instead,  we 
will  consider  iterative  procedures  to  minimize  equation  (3.2). 

So  far  we  assumed  all  parameters  to  be  continuous.  The  coronary  angiogram  is 
digitized,  however.  That  is,  we  are  given  samples  of  the  image  on  a  Cartesian  grid, 
and  the  above  integrals  become  summations.  In  the  following  sections  we  will  look 
at  algorithms  for  finding  the  maximum  likelihood  estimates  of  the  parameters  of 
each  model.  We  will  consider  the  one-dimensional  model  first,  because  it  is  simpler 
and  leads  into  a  better  understanding  of  the  two-dimensional  case. 

3.1  One-Dimensional  Model 

The  one-dimensional  model  of  a  densitometric  profile  along  a  line  crossing  the  coro¬ 
nary  artery  is  described  by  the  discrete  version  of  equation  (2.22) 

f\t]  =  v[f)  *  g[t\  +  p'\t\  +  u/[f]  (3.6) 

The  vessel  component  is  described  by  the  location  of  the  center  c,  the  radius  r,  and 
the  parameter  a.  The  background  component  is  characterized  by  the  coefficients 
of  the  polynomial  dt,  the  Gaussian  point  spread  function  is  characterized  by  the 
standard  deviation  s,  and  the  noise  by  the  standard  deviation  a. 

The  choice  of  the  degree  of  the  polynomial  that  models  the  background  depends 
on  the  width  of  the  vessel  relative  to  the  width  of  the  profile  we  are  analyzing  and 


the  variation  in  the  background.  We  want  the  models  of  the  background  and  vessel 
to  be  independent,  that  is.  the  polynomial  must  fit  the  background  component  not 
the  vessel  component.  Thus  we  must  choose  a  polynomial  that  fits  the  background 
with  as  low  a  degree  as  possible.  Reasonable  values  for  the  degree  of  the  polynomial 
are  three  and  five  (odd  degrees  seem  to  have  better  behavior  than  even  degrees). 

We  can  find  the  maximum  likelihood  estimates  of  the  parameters  character¬ 
izing  our  model  by  maximizing  the  discrete  version  of  the  likelihood  function  of 
equation  (3.2),  or  equivalently  minimizing 

A'(a)  =  £{/[<; -/oi<;aj}2  (3.7) 

4=0 

where 

/o[f;aj  =  vj<!  *  g\t\  -  p"t\  (3.8) 

over  all  parameter  vectors  a. 

Minimizing  A'  over  the  parameters  c,  r,  a,  s ,  and  the  coefficients  of  the  polyno¬ 
mial  is  the  same  as  fitting  a  polynomial  to  all  possible  functions 

fit)  -  i it)  *  g[t\  (3.9) 

and  selecting  the  one  that  results  in  the  minimum  mean  square  error.  This  means 
that  we  must  search  over  all  possible  values  of  the  unknown  parameters.  However, 
taking  the  structure  of  the  problem  into  consideration,  we  can  proceed  with  the 
search  in  a  systematic  way. 

We  assume  that  the  variance  of  the  Gaussian  point  spread  function  can  be 
estimated  independently  of  the  other  parameters.  The  above  minimization  is  done 
as  follows.  First  we  find  crude  estimates  of  the  parameters  c.  r,  and  a.  The 
estimates  of  the  previous  profile  should  be  adequate  because  of  the  spatial  continuity 
of  the  vessel.  At  the  first  profile  we  can  use  some  heuristic  technique  to  determine 
the  estimates.  Once  we  have  the  initial  estimates,  we  need  only  search  in  the 
neighborhood  of  these  estimates.  To  reduce  the  search  even  further,  we  decouple 


the  search  over  c,  and  the  parameters  r  and  a.  This  is  possible,  because  of  the  bell 
shaped  vessel  And  its  model.  It  should  be  obvious  that  with  a  reasonable  estimate 
of  r  and  a,  optimization  over  c  will  give  a  value  pretty  close  to  the  optimal.  We 
can  then  optimize  over  r  and  a.  A  possible  strategy  for  the  minimization  is  the 
following: 

1.  Get  initial  estimates  for  c,  r,  and  a.  Pick  search  increment. 

2.  Update  c  estimate  by  searching  over  a  neighborhood  of  c. 

3.  Update  r  and  a  estimates  by  searching  over  their  neighborhoods. 

4.  Reduce  search  increment. 

5.  Update  c  estimate  by  searching  over  a  neighborhood  of  c. 

6.  Update  r  and  a  estimates  by  searching  over  their  neighborhoods. 

7.  Repeat  last  2  steps  if  necessary. 

It  also  follows  from  the  shape  of  the  vessel  that  the  likelihood  function  is  locally 
convex.  As  we  move  away  from  the  optimum  value  of  one  of  the  parameters  (in 
either  direction)  the  value  of  the  error  should  keep  increasing.  Thus  we  need  only 
search  in  the  direction  of  decreasing  error. 

For  each  value  of  c,  r,  anu  a  we  have  to  evaluate  the  function  in  equation  (3.9) 
and  fit  a  polynomial  to  this  function.  The  main  computation  in  each  function 
evaluation  is  the  convolution  with  the  Gaussian.  For  the  convolution  the  number 
of  computations  is  of  the  order  of  X  log2  N.  The  polynomial  fit  is  of  the  order  of 
QX  computations,  where  Q  is  the  degree  of  the  polynomial.  We  will  refer  to  the 
above  operations  as  one  iteration.  Thus  the  number  of  computations  per  iteration 
is  proportional  to 

computations/iteration  ~  Xlog^N  QN  (3.10) 

Exact  knowledge  of  the  covariance  matrix  of  the  Gaussian  point  spread  function 
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is  not  crucial  in  the  optimization  procedure.  In  fact,  if  we  assumed  that  there  is  no 
blurring.  we'WOuld  still  get  good  estimates  of  the  parameters. 

So  far  we  considered  only  one  profile  of  a  line  crossing  the  vessel.  In  order  to 
apply  the  algorithm  to  a  segment  of  the  vessel,  we  need  a  starting  point  and  the 
approximate  direction  of  the  vessel.  We  assume  that  they  can  be  provided  by  the 
angiographer.  We  can  then  apply  the  algorithm  to  the  profiles  perpendicular  to  the 
given  direction.  The  output  of  the  algorithm  consists  of  the  sequences  of  estimated 
parameters  r  q  ,  a  q\,  and  c  q  .  From  c  q  and  the  direction  of  the  vessel  we  get  the 
axis  coordinates  cl}q  and  ciiq'i. 

However,  the  objective  is  to  estimate  the  dimensions  of  the  cross-sections  that 
are  perpendicular  to  the  vessel  axis.  Using  the  initial  estimate  of  the  axis  coordi¬ 
nates.  we  can  obtain  the  perpendicular  profiles  at  each  point  and  get  new  parameter 
estimates.  This  procedure  can  be  repeated  iteratively.  Usually  one  additional  it¬ 
eration  is  enough.  At  each  point  along  the  vessel  the  perpendicular  direction  is 
computed  by  fitting  a  cubic  spline  to  the  estimated  c  q  sequence,  and  computing 
its  slope.  The  density  of  the  linear  profile  is  obtained  by  simple  linear  interpolation. 

Thus  far  the  one-dimensional  algorithm  does  not  use  the  spatial  continuity  of 
the  vessel,  except  for  determining  initial  estimates  for  the  parameters  of  each  profile 
from  the  previous  one.  In  order  to  have  a  fair  comparison  of  the  one-dimensional 
approach  with  the  two-dimensional  approach,  we  will  consider  a  simple  heuristic  way 
to  use  the  continuity  of  vessel.  This  is  to  smooth  the  vessel  parameter  sequences 
after  they  are  estimated,  by  fitting  splines  to  them.  This  approach  will  be  useful 
in  chapter  5,  where  we  will  compare  the  one-dimensional  approach  to  the  two- 


3.2  Two-Dimensional  Model 


We  consider  a  small  rectangular  section  of  the  angiogram  (N  x  A/.  where  usually 
A /  =  .V  -  50-80  pixels).  We  assume  that  there  is  a  single  isolated  vessel  in  this 
section  of  the  angiogram.  We  assume  that  the  angiographer  can  pick  such  a  section. 
Physical  considerations  dictate  that  the  vessel  must  cross  the  section  borders  in  two 
places  at  most.  It  may  terminate  (total  obstruction)  before  it  reaches  another  part  of 
the  border.  For  simplicity  we  assume  that  the  vessel  does  not  cross  the  two  vertical 
sides  of  the  section.  This  is  not  a  limiting  assumption,  since  the  angiographer  can 
pick  the  section  so  that  it  is  true. 

The  two-dimensional  model  of  the  rectangular  section  is  described  by  the  discrete 
version  of  equation  (2.18) 

f  x.y  =  vjx,  yi  *  g{x,y\  4-  p  |x,y]  -  w{x,y]  (3.11) 

The  vessel  component  is  described  by  the  sequence  of  axis  coordinates  Ci\q\  and  e2'<7  ', 
and  the  sequences  of  the  corresponding  radii  r\q\  and  alphas  algj.  The  background 
component  is  characterized  by  the  coefficients  of  the  polynomial  dX].  Finally,  the 
Gaussian  point  spread  function  is  characterized  by  the  covariance  matrix  5,  and 
the  noise  by  the  variance  o2 . 

We  assume  that  the  variation  in  the  direction  of  the  axis  of  the  vessel  is  limited 
so  that  the  vessel  crosses  each  horizontal  line  at  most  once.  This  is  a  reasonable 
assumption,  given  the  fact  that  we  look  at  a  short  segment  of  the  vessel.  In  this  case 
we  can  choose  c2  q  =  q.  Thus,  the  total  number  of  parameters  describing  the  vessel 
is  3 iV,  where  N  is  the  vertical  dimension  of  the  image  section.  However,  since  we 
have  assumed  that  these  sequences  are  samples  of  continuous  cubic  splines,  they  can 
also  be  described  by  the  spline  parameters.  The  location  of  the  nodes  determines 
how  smooth  the  splines  are.  Thus  we  can  choose  the  spline  parameters  to  reflect 
our  a  priori  knowledge  of  the  vessel  smoothness.  Of  course  the  choice  also  depends 
on  the  magnification  of  the  x-ray. 
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The  choice  of  the  degree  of  the  two-dimensional  polynomial  that  models  the 
background  “de'pends  on  the  size  of  the  vessel  relative  to  the  image  section  we  are 
analyzing  and  the  variation  in  the  background.  As  we  argued  in  the  one-dimensional 
case,  we  want  the  models  of  the  background  and  vessel  to  be  independent,  and  thus 
we  must  choose  a  polynomial  that  fits  the  background  with  as  low  a  degree  as 
possible.  Reasonable  values  for  the  degree  of  the  polynomial  are  three  and  five  (odd 
degrees  seem  to  have  better  behavior  than  even  degrees). 

We  can  find  the  maximum  likelihood  estimates  of  the  parameters  character¬ 
izing  our  model  by  maximizing  the  discrete  version  of  the  likelihood  function  of 
equation  (3.2),  or  equivalently  minimizing 


where 


A'(a)  =  £{/[*.  vl  -  /o[x, y; a]}2 

z=0  y=0 


/0[X,  y; a]  =  u[z,  y]  *  yix,  y]  +  p'\x,  yj 


(3.12) 


(3.13) 


over  all  parameter  vectors  a. 


Minimizing  A'  over  the  model  parameters  is  the  same  as  fitting  a  polynomial  to 
all  possible  functions 

f[x,y]  -  v{x,y\*  g\x,y\  (3.14) 

and  selecting  the  one  that  results  in  the  minimum  mean  square  error.  This  means 
that  we  must  do  an  exhaustive  search  through  all  possible  sequences  cj(y],  r[y]  and 
a  y  .  This  is  a  formidable  task.  However,  as  we  did  in  the  one-dimensional  case, 
we  can  take  the  structure  of  the  problem  into  consideration  and  proceed  with  the 
search  in  a  systematic  way. 

We  assume  that  the  covariance  matrix  of  the  Gaussian  point  spread  function 
can  be  estimated  independently  of  the  other  parameters.  The  above  minimization 
is  done  a s  follows.  First  we  determine  a  crude  estimate  of  the  sequences  Ci[yj,  r[q], 
and  and  then  we  consider  perturbations  around  these  estimates.  For  this  we 
can  use  the  estimates  obtained  by  some  heuristic  technique,  or  the  one-dimensional 


r~W  rm 


approach  presented  in  the  previous  section.  We  can  also  obtain  initial  estimates  of 
the  vessel  directly  by  a  technique  described  in  section  3.3. 

To  simplify  the  search,  we  decouple  the  search  over  the  sequences  cj  jgj,  r  q  and 
o [q\.  As  we  argued  in  the  one-dimensional  case,  this  is  possible  because  of  the  bell 
shaped  vessel  and  its  model.  It  also  follows  from  the  shape  of  the  vessel  that  the 
likelihood  function  is  convex  for  perturbations  of  a  sequence  around  a  point.  One 
possible  strategy  for  the  minimization  is  the  following: 

1.  Get  initial  estimate  of  the  sequences  r[<?j,  and  a[g). 

2.  Optimize  over  a  set  of  rough  perturbations  around  the  current  r\q j  estimates. 

3.  Optimize  over  a  set  of  rough  perturbations  around  the  current  cj[g]  estimates. 

4.  Optimize  over  a  set  of  fine  perturbations  around  the  current  r[g]  estimates. 

5.  Optimize  over  a  set  of  fine  perturbations  around  the  current  Ci [<7 j  estimates. 

6.  Consider  tradeoff  in  perturbations  around  the  r[<?]  estimates  and  alg]  esti¬ 
mates. 

7.  Repeat  last  3  steps  if  necessary. 

At  each  step,  we  consider  a  set  of  perturbations  around  a  parameter  sequence. 
For  example  consider  the  sequence  r[qrj  shown  in  figure  3.1.  The  solid  line  shows  the 
original  sequence,  and  the  dotted  line  shows  a  typical  perturbation  of  step  4.  The 
perturbations  are  made  at  K  evenly  spaced  points,  shown  by  the  vertical  dashed 
lines  in  the  figure.  As  can  be  seen  in  the  picture,  the  perturbations  are  local. 
Our  ability  to  do  this  follows  from  the  flexibility  of  the  splines.  As  we  discussed  in 
chapter  2,  the  splines  can  be  changed  locally  without  affecting  their  global  behavior. 
On  the  other  hand  the  properties  of  the  splines  guarantee  that  these  changes  are 
smooth.  Thus  the  properties  of  the  splines  insure  that  spatial  continuity  constrains 
the  vessel  parameters  only  locally. 

Another  property  of  the  splines  that  is  important  here  is  the  fact  that  their 
flexibility  can  be  adjusted  by  changing  the  spacing  of  the  nodes.  Figure  3.2  shows 


Figure  3.1:  Example  of  Fine  Perturbation 

an  example  of  a  perturbation  at  an  earlier  stage  in  the  optimization  procedure 
(step  2).  Here  the  spacing  of  the  nodes  is  wider  and  the  spline  less  flexible.  As  a 
result  the  spacing  of  the  points  at  which  the  perturbations  are  made  need  not  be 
as  dense. 

For  each  perturbation  we  have  to  evaluate  the  function  in  equation  (3.14)  and  fit 
a  polynomial  to  this  function.  Each  function  evaluation  involves  the  computation 
of  the  vessel  component,  the  convolution  with  a  Gaussian,  and  a  subtraction.  The 
number  of  computations  for  the  vessel  evaluation  is  of  the  order  of  N2  (from  now  on 
we  set  M  =  N).  For  the  convolution  the  number  of  computations  is  of  the  order  of 
N2  log2  N.1  Finally,  the  polynomial  fit  is  of  the  order  of  \Q2N2  computations,  where 
Q  is  the  degree  of  the  polynomial.  We  will  refer  to  the  above  cycle  of  computations 
as  one  iteration. 

We  observe  that  the  computations  are  dominated  by  the  convolution  with  the 

‘Actually  the  number  of  computations  is  of  the  order  of  N2  log2  jVj  ,  where  N\  >  ,V  ■+•  NtJ  and 
.V„  x  N„  is  the  support  of  the  Gaussian  point  spread  function. 
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Figure  3.2:  Example  of  Crude  Perturbation 

Gaussian  point  spread  function  and  the  polynomial  fitting.  Considering  only  these 
two  operations,  the  number  of  computations  per  iteration  is  proportional  to 

computations/iteration  ~  N2log2N  -r  - Q2N 2  (3.15) 

2 

Exact  knowledge  of  the  covariance  matrix  of  the  Gaussian  point  spread  function 
is  not  crucial  in  the  optimization  procedure.  In  fact  we  can  assume  that  there  is 
no  blurring  until  the  very  last  optimization  stage,  without  significantly  affecting 
the  procedure.  Thus  we  can  save  a  considerable  amount  of  computation,  since  the 
number  of  operations  per  iteration  becomes  proportional  to  Q2N 2 . 

3.2.1  Comparison  With  One-Dimensional  Approach 

In  terms  of  computations,  the  advantage  of  the  one-dimensional  approach  is  that 
it  optimizes  the  cross-section  parameters  of  the  vessel  at  each  point  independently 
of  all  others.  Also,  it  only  uses  background  information  of  one  linear  profile.  The 
two-dimensional  optimization  is  also  done  locally,  but  more  than  one  cross-section  is 


involved  at  each  point,  as  the  algorithm  takes  the  spatial  continuity  of  the  vessel  into 
account.  Also. 'it  uses  all  of  the  background,  and  assumes  blurring  in  all  directions. 
Thus  we  expect  the  two-dimensional  approach  to  require  much  more  computation. 
This  is  indeed  the  case  as  demonstrated  by  the  following  crude  comparison. 

The  number  of  computations  at  each  step  of  the  two-dimensional  optimization 
procedure  is  proportional  to  K ,  the  number  of  perturbation  points,  and  the  num¬ 
ber  of  computations  per  iteration  of  equation  (3.15).  Thus  the  total  number  of 
computations  is  proportional  to 

computations  for  2D  case  ~  K(\2log2N  -  ~Q2\2)  (3.16) 

The  number  of  computations  at  each  step  of  the  one-dimensional  optimization 
procedure  is  proportional  to  the  number  of  computations  per  iteration  of  equa¬ 
tion  (3.10).  But  this  is  repeated  for  ,V  profiles  along  the  vessel.  Thus  the  total 
number  of  computations  is  proportional  to 

computations  for  ID  case  ~  .X2log2N  -  QN2  (3.17) 

Assuming  that  everything  else  is  the  same,  we  can  see  that  the  two-dimensional 
approach  is  much  more  time  consuming  than  the  one-dimensional  approach.  In 
fact,  the  two-dimensional  approach  requires  more  computation  than  this  simple 
comparison  indicates,  mainly  because  more  iterations  are  needed.  Finally,  it  should 
be  obvious  that  both  approaches  involve  much  more  computation  than  the  slope 
methods,  mainly  because  of  the  search  involved  in  the  optimization  procedure. 

3.3  Initial  Vessel  Estimate 

The  modeling  approach  presented  in  this  thesis  suggests  a  way  to  estimate  the 
background  from  the  angiogram  without  estimating  the  vessel  parameters.  At  the 
same  time  we  can  get  a  rough  estimate  of  the  vessel  in  the  section. 


Consider  a  rectangular  section  of  the  coronary  angiogram,  containing  one  or 
more  vesselsTAssume  that  the  background  is  slowly  varying  compared  to  the  vessel 
profile,  and  thus  can  be  modeled  by  a  low  order  two-dimensional  polynomial.  This 
assumption  is  reasonable,  as  we  discussed  in  the  previous  chapter,  if  the  section  of 
the  angiogram  to  be  analyzed  is  chosen  properly.  Also  assume  that  the  vessel  profile 
is  positive,  which  is  always  true.  Now  consider  the  following  algorithm. 

Algorithm:  Background  Estimation 

1.  Fit  a  low  order  polynomial  to  the  section,  and  compute  the  standard  deviation. 

2.  Remove  the  points  for  which  the  error  is  greater  than  a  given  multiple  of  the 
standard  deviation. 

3.  Fit  a  low  order  polynomial  to  the  remaining  points  of  the  section,  and  compute 
the  standard  deviation  35j. 

4.  Quit  if  there  is  no  significant  reduction  in  the  standard  deviation.  Otherwise 
go  back  to  step  2. 

This  algorithm  has  been  observed  to  converge  after  a  few  iterations.  The  polyno¬ 
mial  thus  obtained  is  a  good  estimate  of  the  background.  The  set  of  points  that 
are  removed  in  the  process  give  a  good  initial  estimate  of  the  vessel  after  we  use 
some  simple  neighborhood  operations  to  remove  isolated  points  due  to  noise  in  the 
background. 

An  example  of  this  procedure  is  shown  in  figure  3.3.  The  original  coronary 
angiogram  section  is  shown  in  figure  3.3(a).  After  3  iterations  of  algorithm  1,  the 
background  estimate  is  shown  in  figure  3.3(b).  The  degree  of  the  background  poly¬ 
nomial  used  in  this  example  is  3,  and  the  threshold  used  was  1.5  times  the  standard 
deviation.  The  set  of  points  above  the  threshold  are  shown  in  figure  3.3(c).  Fig¬ 
ure  3.3(d)  shows  the  points  after  removal  of  isolated  points  by  simple  neighborhood 
operations. 


Chapter  4 


Comparison  of  One-Dimensional 
Model  With  Slope  Methods 

In  this  chapter  we  examine  the  performance  of  the  one-dimensional  model  and 
compare  it  to  that  of  the  slope  methods,  on  which  most  current  boundary  detection 
methods  are  based.  In  our  analysis  we  will  consider  each  densitometric  profile  in¬ 
dependently  of  the  previous  profiles,  i.e..  the  spatial  continuity  of  the  vessel  or  the 
background  (in  the  direction  of  the  vessel)  will  not  be  taken  into  account.  Our  objec¬ 
tive  is  to  isolate  the  problem  of  accurately  estimating  the  cross-sectional  dimensions 
from  one  profile,  from  that  of  exploiting  the  spatial  continuity  of  the  vessel.  The 
importance  of  vessel  continuity  will  be  examined  in  the  following  chapter. 

The  slope  method  we  implemented  in  our  examples  is  the  following:  at  each  point 
of  the  profile  we  fitted  a  second  degree  polynomial  to  a  set  of  7  adjacent  points. 
The  first  and  second  derivative  of  the  profile  were  estimated  as  the  corresponding 
derivatives  of  the  polynomial.  The  inflection  points  were  identified  as  the  points 
with  maximum  and  minimum  first  derivative,  and  the  knee  points  as  the  maxima  of 
the  second  derivative.  Variations  of  the  above  approach,  as  well  as  other  approaches 
of  estimating  the  profile  derivatives,  were  tried  with  comparable  or  worse  results. 

The  vessel  diameter  was  obtained  as  the  difference  between  the  inflection  or  the 


knee  points.  The  cross-sectional  .  5a  was  obtained  as  the  integral  of  the  profile 
between  the-boundary  points  after  the  background  had  been  estimated  and  sub¬ 
tracted.  The  background  was  estimated  by  computing  the  linear  regression  line 
through  the  boundary  points. 

We  tested  our  algorithm  on  computer  generated  data  (synthetic  data),  on  x-rays 
of  contrast-medium-filled  cylindrical  phantoms,  and  on  real  coronary  angiograms. 
The  first  two  cases  are  useful,  because  the  results  can  be  compared  to  the  true 
vessel  dimensions.  The  last  case  is  useful  as  an  illustration  of  application  of  the 
algorithms  to  actual  coronary  angiograms.  The  true  dimensions  are  not  known,  so 
any  conclusions  we  draw  will  be  only  qualitative. 

In  the  phantom  and  coronary  artery  cases  a  sequence  of  images  were  recorded 
on  35mm  film  (cine-angiogram)  and  a  single  frame  was  selected  by  the  angiogra- 
pher.  The  analog  image  was  optically  magnified  and  digitized  by  a  vidicon  cam¬ 
era/digitizer  interfaced  to  a  Vax  11/780  computer.  The  grey  level  of  each  pixel  was 
digitized  to  256  levels.  The  scans  were  performed  four  times  and  averaged.  The 
phantom  x-rays  were  the  same  ones  used  by  Spears  et  a/,  in  [32] ;  the  resolution 
was  fixed  at  approximately  15  microns/pixel.  The  resolution  for  the  real  coronary 
angiograms  varied  from  image  to  image.  In  the  synthetic  data  case  the  grey  level 
of  each  pixel  was  also  digitized  to  256  levels  to  be  consistent  with  the  x-ray  data. 

4.1  Synthetic  Data 

First  we  examine  the  performance  of  the  model  on  synthetic  data.  The  data  were 
generated  assuming  the  model  is  valid.  A  typical  example  of  a  synthetic  densito- 
metric  profile  is  shown  in  figure  4.1.  This  profile  is  80  pixels  wide.  The  background 
is  a  fifth  degree  polynomial,  obtained  from  a  real  angiogram  profile.  The  vessel 
parameter  values  are  c  =  40  pixels,  r  =  10  pixels,  a  =  2,  and  s  =  1  pixel.  These 
parameters  were  chosen  to  correspond  to  vessel  profiles  of  actual  angiograms.  Fig- 


Figure  4.1:  Synthetic  Densitometric  Profile 

ures  4.1(a)  and  (b)  show  the  profile  before  and  after  the  noise  was  added.  The 
standard  deviation  of  the  noise  is  a  —  2,  which  corresponds  to  signal  to  noise  ratio 
of  19.6  dB. 

We  conducted  a  number  of  experiments  to  study  the  properties  of  the  modeling 
approach  and  to  compare  them  to  the  properties  of  the  slope  method.  In  particu¬ 
lar,  we  investigated  the  behavior  of  the  two  approaches  for  different  levels  of  noise, 
different  vessel  dimensions,  and  different  image  backgrounds.  In  each  of  the  ex¬ 
periments  we  describe  below,  all  the  parameters  (background,  c,  r,  a,  s,  and  noise 
variance)  were  fixed  except  one.  We  considered  80  noise  samples,  and  estimated 
the  vessel  dimensions  for  each  sample  using  the  modeling  approach  and  the  slope 
methods.  The  radius  and  cross-sectional  area  estimates  were  averaged,  and  the 
standard  deviation  of  the  estimates  was  computed.  The  standard  deviation  of  the 
estimates  was  normalized,  that  is,  it  was  divided  by  the  true  value  of  the  parameter 
being  estimated  (radius  or  area). 

The  degree  of  the  polynomials  used  in  the  modeling  approach  was  5,  the  same  as 
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that  of  the  background  polynomials  of  the  synthetic  data.  In  all  of  the  experiments 
the  model  parameters  were  chosen  to  correspond  to  profiles  of  actual  angiograms. 


4.1.1  Effect  of  Noise  Level 

For  this  experiment  all  the  parameters  were  fixed  except  the  noise  level.  The  syn¬ 
thetic  densitometric  profile  we  used  was  that  of  figure  4.1(a).  The  vessel  parameters 
were  c  =  40  pixels,  r  =  10  pixels,  a  =  2.  and  s  =  1  pixel.  To  this  profile  we  added 
7  different  levels  of  noise.  For  each  noise  level  we  considered  80  noise  samples  and 
estimated  the  vessel  dimensions  for  each  sample.  The  signal  to  noise  ratio  (SNR) 
is  defined  as  follows: 


SNR  =  10 


log 


10 


signal  variance\ 
noise  variance  ) 


(4.1) 


The  average  of  the  radius  estimates  as  a  function  of  signal  to  noise  ratio  is  shown 
in  figure  4.2(a)  and  the  normalized  standard  deviation  is  shown  in  figure  4.2(b). 
The  solid  fine  represents  the  model  estimates,  the  dot-dashed  line  represents  the 
corresponding  results  when  the  inflection  points  are  found,  and  the  dashed  line 
represents  the  results  for  the  knee  points.  The  true  radius  is  shown  in  dotted  line. 
We  observe  that  the  inflection  point  method  underestimates  the  vessel  radius,  while 
the  knee  point  method  overestimates  it.  This  behavior  is  characteristic  of  such 
algorithms,  and  is  consistent  with  previous  results  [321.  Note  that  for  the  signal 
to  noise  ratios  considered,  the  average  of  the  model  estimates  is  roughly  constant. 
The  same  is  true  for  the  inflection  point  estimates,  but  the  average  of  the  knee 
point  estimates  increases  for  low  SNR’s.  As  expected  the  standard  deviation  of  the 
estimates  is  significantly  higher  in  the  slope  methods. 

The  average  of  the  area  estimates  as  a  function  of  signal  to  noise  ratio  is  shown 
in  figure  4.3(a)  and  the  standard  deviation  is  shown  in  figure  4.3(b).  The  solid 
line  represents  the  model  estimates,  the  dot-dashed  line  shows  the  inflection  point 
estimates,  and  the  dashed  line  shows  the  knee  point  estimates.  The  dotted  line 
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Figure  4.2:  Radius  Estimation  for  Different  SNR’s 
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Figure  4.3:  Area  Estimation  for  Different  SNR’s 


represents  the  true  area.  Here  we  observe  that  the  knee  point  method  slightly 
underestimates  the  area,  but,  more  importantly,  the  estimate  changes  significantly 
with  the  signal  to  noise  ratio.  For  the  signal  to  noise  ratios  considered,  the  average 
of  the  model  estimates  is  roughly  constant.  The  inflection  point  method  grossly 
underestimates  the  area.  As  a  result,  the  usefulness  of  the  inflection  point  area 
estimates  is  doubtful,  but  we  include  the  results  for  completeness.  The  standard 
deviation  of  the  slope  approach  estimates  is  higher  than  that  of  the  model  estimates. 

4.1.2  Effect  of  Vessel  Radius 

For  this  experiment  all  the  parameters  were  fixed  except  the  vessel  radius.  The 
vessel  parameters  were  c  —  40  pixels,  a  =  2,  and  s  =  1  pixel.  The  standard  deviation 
of  the  noise  was  fixed  at  a  =  2  (this  corresponds  to  signal  to  noise  ratio  of  19.6  dB 
when  r  =  10  pixels).  The  values  of  the  radius  were  r  =  3,4,5,6,7,8,10,12,15 
pixels.  For  each  radius  we  considered  80  noise  samples  and  estimated  the  vessel 
dimensions  for  each  sample. 

The  average  of  the  radius  estimates  as  a  function  of  radius  is  shown  in  fig¬ 
ure  4.4(a)  and  the  normalized  standard  deviation  is  shown  in  figure  4.4(b).  The 
solid  line  represents  the  model  estimates,  the  dot-dashed  line  represents  the  inflec¬ 
tion  point  estimates,  and  the  dashed  line  represents  the  knee  point  estimates.  The 
true  radii  lie  on  the  diagonal  (dotted  line).  Observe  that  the  inflection  point  method 
underestimates  the  vessel  radius,  and  the  knee  point  method  overestimates  it.  The 
problem  for  the  slope  methods,  however,  is  that  the  estimates  lie  on  a  line  that  does 
not  pass  through  the  origin.  Moreover,  the  line  curves  up  near  the  origin.  This 
means  that  the  relative  radius  estimates  are  biased,  and  thus  can  lead  to  significant 
errors  in  the  determination  of  vessel  obstruction.  We  also  observe  that  the  standard 
deviation  of  ’he  slope  approach  estimates  is  higher. 

The  average  of  the  area  estimates  as  a  function  of  area  is  shown  in  figure  4.5(a) 
and  the  normalized  standard  deviation  is  shown  in  figure  4.5(b).  Both  axes  of 
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Figure  4.5:  Area  Estimation  for  Different  Radii 


figure  4.5(a)  and  the  x-axis  of  figure  4.5(b)  show  the  square  root  of  the  area,  so 
that  the  data  points  are  more  evenly  spaced.  The  solid  line  represents  the  model 
estimates,  the  dot-dashed  line  represents  the  inflection  point  estimates,  and  the 
dashed  line  represents  the  knee  point  estimates.  The  true  areas  lie  on  the  diagonal 
(dotted  line).  Here  the  knee  point  estimates  of  the  area  lie  on  a  line  very  close  and 
roughly  parallel  to  the  diagonal.  Again  the  line  bends  up  near  the  origin.  These 
observations  are  much  more  pronounced  in  the  inflection  point  results.  We  also 
observe  that  the  standard  deviation  of  the  slope  approach  estimates  is  higher.  The 
average  of  the  model  estimates  is  on  the  diagonal. 

4.1.3  Effect  of  Parameter  Alpha 

For  this  experiment  all  the  parameters  were  fixed  except  a.  The  vessel  parameters 
were  c  =  40  pixels,  r  =  8  pixels,  and  s  =  1  pixel.  The  standard  deviation  of  the 
noise  was  fixed  at  a  —  2  (this  corresponds  to  signal  to  noise  ratio  of  17.2  dB  when 
a  =  l).  The  parameter  alpha  ranged  from  a  =  0.5  to  o  =  2.5  in  steps  of  0.5.  For 
each  alpha  we  considered  80  noise  samples  and  estimated  the  vessel  dimensions  for 
each  sample. 

The  average  of  the  radius  estimates  as  a  function  of  alpha  is  shown  in  fig¬ 
ure  4.6(a)  and  the  normalized  standard  deviation  is  shown  in  figure  4.6(b).  The 
solid  line  represents  the  model  estimates,  the  dot-dashed  line  represents  the  inflec¬ 
tion  point  estimates,  and  the  dashed  line  represents  the  knee  point  estimates.  The 
true  radius  is  shown  in  dotted  line.  Again  we  observe  that  the  inflection  point 
method  underestimates  the  vessel  radius,  and  the  knee  point  method  overestimates 
it.  Note  that  the  knee  point  estimates  increase  for  small  alphas.  The  standard  de¬ 
viation  of  the  slope  approach  estimates  is  higher  than  that  of  the  model  estimates. 

The  average  of  the  area  estimates  as  a  function  of  alpha  is  shown  in  figure  4.7(a) 
and  the  normalized  standard  deviation  is  shown  in  figure  4.7(b).  The  solid  line 
represents  the  model  estimates,  the  dot-dashed  line  represents  the  inflection  point 


Radius  [pixels] 


- 1  - - 1' 


inflection  points 


knee  points 


3 

I* 

s* 


Alpha  [ ] 


*  s' 
V 

•  V 


(a)  Average  of  Radius  Estimates 


Standard  Deviation  [normalized] 


0.15- 


Alpha  [  ] 


(b)  Normalized  Standard  Deviation  of  Radius  Estimates 


Figure  4.6:  Radius  Estimation  for  Different  Alphas 


Are*  [plxels-2] 


(a)  Average  of  Area  Estimates 


Standard  Deviation  [normalized] 


Inflection  points 


knee  points 


Area  [pixels-2] 


k  '  U 

f  - 


Area  [pixels-2] 


(b)  Normalized  Standard  Deviation  of  Area  Estimates 


Figure  4.7:  Area  Estimation  for  Different  Alphas 


estimates,  and  the  dashed  line  represents  the  knee  point  estimates.  The  true  areas 
lie  on  the  diagonal  (dotted  line).  Here  we  observe  a  significant  bias  for  the  knee 
point  method.  The  estimates  lie  on  a  line  that  curves  up  as  we  approach  the  origin. 
Also  the  standard  deviation  of  the  slope  approach  estimates  is  higher  than  that  of 
the  model  estimates. 

4.1.4  Effect  of  Background 

For  this  experiment  all  the  parameters  were  fixed  except  the  background.  The  vessel 
parameters  were  c  -  40  pixels,  r  =  10  pixels,  a  -  1,  and  s  -  1  pixel.  The  standard 
deviation  of  the  noise  was  fixed  at  a  -  2  (this  corresponds  to  signal  to  noise  ratio 
of  10.9  dB  for  the  flat  background).  The  different  backgrounds  used  are  shown 
in  figure  4.8.  They  were  obtained  from  actual  angiogram  profiles.  For  each  back¬ 
ground  we  considered  80  noise  samples  and  estimated  the  vessel  dimensions  for  each 
sample.  The  average  of  the  radius  estimates  for  the  different  backgrounds  is  shown 
in  figure  4.9(a)  and  the  normalized  standard  deviation  is  shown  in  figure  4  9(b) 
The  solid  line  represents  the  model  estimates,  the  dot-dashed  line  represents  the 
corresponding  results  when  the  inflection  points  are  found,  and  the  dashed  line 
represents  the  results  for  the  knee  points.  The  true  radius  lies  on  the  dotted  line 
The  corresponding  results  for  the  area  estimates  are  shown  in  figures  4  101  a)  and 
4.10(b).  We  observe  that  the  model  estimates  are  not  very  sensitive  to  variations  m 
the  background.  The  same  is  true  for  the  radius  estimates  of  the  slope  approaches 
but  not  for  the  area  estimates.  This  is  expected,  because  the  background  is  smooth 
and  thus  cannot  have  a  significant  effect  on  the  derivatives  of  the  profile,  but  it  has 
an  effect  on  the  area  estimates  which  assume  that  the  background  is  linear  Con¬ 
sistent  with  the  previous  experiments,  the  standard  deviation  of  the  slope  approach 
estimates  is  higher  than  that  of  the  model  estimates 
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Figure  4.9:  Radius  Estimation  for  Different  Backgrounds 


Area  [plxals-2] 


5 


V 

v 


.V 

t 


$ 


600  r 


400  1 


H - \ 


200- 


t" 


1 _ L 


(a)  Average  of  Area  Estimates 


Standard  Deviation  [normalized] 

0.2: - - - I— 


i - r 


0.15- 


0.1  - 


0.05' 


V. 


'■r 

Y 


k 

v _ L 

1 - '  ^ 


I - - 


model 


Inflection  points 
knee  points 
true 


H 


0  1  2  3  4  5  6  7 


Background 


H 


Background 


(b)  Normalized  Standard  Deviation  of  Area  Estimates 


Figure  4.10:  Area  Estimation  for  Different  Backgrounds 
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4.1.5  Remarks 


Summarizing  the  results  of  the  synthetic  data  experiments,  we  note  that  in  all 
cases  the  performance  of  the  modeling  approach  is  superior  to  that  of  the  slope 
methods.  It  presents  an  improvement  in  both  the  average  of  the  estimates  and  in 
their  standard  deviation.  The  bias  in  the  slope  methods  can  lead  to  significant  errors 
in  the  determination  of  the  relative  dimensions  at  different  points  of  the  vessel. 

Our  results  also  indicate  that  the  knee  point  estimates  are  less  reliable  than 
the  inflection  point  estimates.  This  should  be  expected,  because  the  knee  point 
estimates  use  the  second  derivatives  of  the  profiles,  whereas  the  inflection  point 
estimates  use  the  first  derivatives,  which  are  less  sensitive  to  noise. 


4.2  Vessel  Phantoms 

Now  we  examine  the  performance  of  the  algorithms  on  the  vessel  phantom  x-rays. 
The  vessel  phantoms  are  cylindrical  holes  with  14  different  radii,  drilled  in  a  plastic 
block,  and  filled  with  angiographic  contrast  material.  The  radiographic  properties  of 
the  plastic  block  are  very  similar  to  those  of  water,  and  thus  the  human  tissue.  The 
block  is  parallelepiped,  so  the  background  of  the  phantom  x-rays  is  flat.  All  x-ravs 
were  obtained  under  the  same  set  of  radiographic  conditions.  The  vessel  phantoms 
are  shown  in  figure  4.11.  Figure  4.12  shows  the  perpendicular  cross-section  of  a 
vessel  phantom  and  the  profile  of  the  perfect  projection  at  that  cross-section.  Since 
the  phantoms  have  constant  cross-sections,  this  experiment  is  similar  to  the  radius 
experiment  performed  in  the  synthetic  data  case.  We  considered  50  consecutive 
profiles  of  each  phantom,  and  estimated  the  vessel  dimensions  for  each  profile.  The 
estimates  were  averaged  and  the  standard  deviation  of  the  estimates  was  computed 
The  degree  of  the  polynomials  used  in  the  modeling  approach  was  3.  As  we  discussed 
in  section  3.1.  a  lower  degree  (3  instead  of  5)  polynomial  to  model  an  approximated 
flat  background  performs  better  in  this  case. 


Figure  4.12:  Vessel  Phantom  Cross-Section  and  Perfect  Projection 


The  average  of  the  radius  estimates  as  a  function  of  true  radius  is  shown  in 
figure  1  13(a)  and  the  normalized  standard  deviation  is  shown  in  figure  4  13(b). 
The  solid  line  represents  the  model  estimates,  the  dot-dashed  line  represents  the 
inflection  point  estimates,  and  the  dashed  line  represents  the  knee  point  estimates 
The  units  of  the  \-axis  are  millimeters  Imrnl  and  of  the  y-axis  pixels  The  exact 
<  orrespondence  is  not  known,  and  is  important  only  if  absolute  measurements  are 
desired  If  we  are  interested  in  relative  measurements,  then  the  estimates  are  satis¬ 
factory  when  thev  lie  on  a  straight  line  pa.ss.ng  through  ’he  origin  Since  the  exact 
correspondence  between  the  radius  in  mm  and  the  radius  m  pixels  is  not  known, 
we  normalized  the  standard  deviation  of  the  estimate^  bv  the  values  on  a  diagonal 
line  (  lose  to  the  average  of  the  estimates 

Observe  that  the  average  of  the  model  estimates  roughly  lies  on  a  straight  line 
through  the  origin  This  is  not  true  for  the  ir  Meet  ion  and  knee  point  estimates 
I  he  i  orrespondtng  lines  do  not  pass  through  ’fie  origin  and  curve  up  near  the 
origin  Observe  also  that  the  inflection  point  method  underest miates  the  vessel 
rad  u*  relative  to  the  modeling  approach  and  the  knee  point  method  overestimates 
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it.  These  results  are  si  nilar  to  those  obtained  in  the  radius  experiment  of  the 
synthetic  data  case.  The  standard  deviation  of  the  model  estimates  is  only  slightly 
lower  than  that  of  the  slope  estimates.  This  is  due  to  the  fact  that  in  this  set  of 
experiments  the  signal  to  noise  ratio  is  very  high  for  phantoms  with  radii  greater 
than  1mm.  Also,  the  fact  that  the  background  is  flat  eliminates  one  potential  source 
of  error  for  the  slope  methods,  which  estimate  the  background  as  the  straight  line 
connecting  the  boundary  points. 

The  average  of  the  area  estimates  as  a  function  of  area  is  shown  in  figure  4.14(a) 
and  the  standard  deviation  is  shown  in  figure  4.14(b).  Similarly  to  the  radius  case, 
the  standard  deviation  of  the  estimates  was  normalized  by  the  values  on  a  diagonal 
line  close  to  the  average  of  the  estimates.  Both  axes  of  figure  4.14(a)  and  the  x-axis 
of  figure  4.14(b)  show  the  square  root  of  the  area,  so  that  the  data  points  are  more 
evenly  spaced.  The  solid  line  represents  the  model  estimates,  the  dot-dashed  line 
represents  the  inflection  point  estimates,  and  the  dashed  line  represents  the  knee 
point  estimates. 

We  observe  that  the  knee  point  method  slightly  underestimates  the  area  relative 
to  the  model  estimates.  As  we  noted  in  the  synthetic  data  experiments,  the  inflection 
point  method  grossly  underestimates  the  area,  and  is  of  limited  usefulness.  Also, 
the  standard  deviation  of  the  model  estimates  is  lower  than  that  of  the  knee  point 
estimates.  However,  there  is  a  major  difference  from  the  synthetic  data  results.  The 
average  of  the  estimates  does  not  lie  on  a  straight  line.  Instead,  it  lies  on  a  line  that 
bends  down  as  the  radius  increases.  This  is  true  for  both  the  modeling  approach 
and  the  slope  approach.  The  reason  for  this  is  that  the  relationship  between  the 
density  of  the  radiographic  image  and  the  thickness  of  the  contrast  agent  is  not 
linear.  As  we  discussed  in  the  introduction,  appropriate  corrections  can  be  made 
to  restore  linearity  30,31,17  .  This  result  indicates  that  the  cross-sectional  area 
estimation  is  more  sensitive  to  the  imaging  process,  and  emphasizes  the  importance 
of  densitometric  checks  and.  if  necessary,  calibration  of  the  imaging  system. 
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Figure  4.14:  Phantom  Experiment:  Area  Estimation 


4.3  Coronary  Angiograms 


The  performance  of  the  algorithms  was  also  tested  on  real  coronary  angiograms. 
As  we  mentioned  earlier,  it  is  very  difficult  to  evaluate  the  performance  of  the 
algorithms  in  this  case,  because  we  do  not  know  the  true  dimensions  of  the  arteries. 
However,  a  simple  inspection  of  the  angiograms  and  their  models  can  give  us  a  good 
indication  of  the  relative  performance  of  the  two  approaches.  Thus  the  comparisons 
in  this  section  will  be  qualitative,  in  contrast  to  the  quantitative  results  of  the 
previous  sections. 

In  the  examples  we  present,  the  approximate  direction  of  the  blood  vessels  is 
vertical.  The  one-dimensional  algorithm  was  applied  to  the  profiles  perpendicular 
to  this  direction.  The  degree  of  the  polynomials  used  in  the  modeling  approach  was 
5.  As  we  discussed  in  section  3.1,  a  second  iteration  of  the  algorithms  could  provide 
the  parameters  of  the  perpendicular  vessel  cross-sections.  Since  both  approaches 
use  the  same  profiles,  we  chose  not  to  carry  out  the  second  step.  This  does  not  affect 
the  comparison  of  the  two  methods,  and  makes  the  presentation  of  the  results  a  lot 
easier. 

A  section  of  a  coronary  angiogram  is  shown  in  figure  4.15(a):  the  resolution 
is  80  x  80  pixels.  The  boundaries  detected  by  the  modeling  approach  are  shown 
in  figure  4.15(b).  Figures  4.15(c)  and  (d)  show  the  inflection  and  knee  points, 
respectively.  The  boundaries  detected  by  the  modeling  approach  and  the  inflection 
points  appear  to  follow  the  shape  of  the  vessel  very  well.  The  knee  point  boundaries 
are  more  noisy,  as  should  be  expected  from  our  analysis  and  previous  results.  Also 
consistent  with  our  previous  results  is  the  fact  that  the  inflection  point  boundaries 
are  inside  the  model  boundaries,  and  the  knee  point  boundaries  are  outside.  As  we 
mentioned  in  section  1.2.  the  inflection  points  are  close  to  what  our  eyes  perceive  as 
the  vessel  edge,  which  is  not  the  same  as  the  actual  vessel  edges.  The  area  estimates 
are  shown  in  figure  4.16.  The  solid  line  represents  the  modeling  approach  estimates, 
the  dot-dashed  line  represents  the  inflection  point  estimates,  and  the  dashed  line 
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Figure  4.18:  Cross-sectional  Area  Detection  (Example  B) 

4.4  Summary 

The  experimental  results  we  presented  in  this  chapter  indicate  that  the  performance 
of  the  modeling  approach  is  better  than  that  of  the  slope  methods.  Our  results 
show  that  the  inflection  points  tend  to  be  inside  the  model  based  boundaries  and 
the  knee  points  tend  to  be  outside.  More  importantly,  this  bias  is  not  consistent, 
thus  resulting  in  errors  in  the  determination  of  the  relative  dimensions  at  different 
points  of  the  vessel.  The  variance  of  the  slope  method  estimates  was  shown  to 
be  higher  than  that  of  the  model  based  estimates.  In  the  case  of  the  knee  point 
estimates,  the  variance  was  especially  high  for  low  signal  to  noise  ratios,  that  is,  for 
the  obstructed  sections  of  the  vessel. 

The  reason  for  the  better  performance  of  the  modeling  approach  is  that  it  uses  all 
•he  information  available  in  the  vessel  profile.  The  slope  approaches  use  information 
•  t, <«allv  (the  slope  at  the  sides,  or  the  second  derivative  near  the  base  points 
1  •tie  profile i  Thus  they  can  be  very  sensitive  to  the  local  shape  of  the  profile. 
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rather  than  its  overall  structure  For  example,  if  there  are  two  local  maxima  of 
the  slope,  it  is  not  clear  which  one  should  be  the  boundary  point.  Similarly,  the 
cross-sectional  area  estimates  depend  only  on  the  boundary  points  to  estimate  the 
background.  The  model  works  well  even  when  the  number  of  samples  across  the 
vessel  profile  is  very  small.  This  is  an  advantage  over  the  slope  methods  in  which 
it  becomes  difficult  to  determine  the  inflection  or  knee  points  accurately  as  the 
number  of  samples  decreases. 

The  advantage  of  the  slope  methods  is  their  computational  simplicity.  However, 
when  one  considers  the  time  and  effort  that  is  required  to  get  an  angiogram  of  the 
coronary  arteries,  it  appears  to  be  reasonable  to  spend  more  time  to  get  better 
estimates  of  the  severity  of  the  obstructions. 


Chapter  5 


Comparison  of  Two-Dimensional 
Model  With  One-Dimensional 
Model 


In  this  chapter  we  examine  the  performance  of  the  two-dimensional  model  and 
compare  it  to  that  of  the  one-dimensional  model  The  two-dimensional  model  in¬ 
corporates  the  spatial  continuity  of  the  vessel  and  the  background  and  is  thus 
expected  to  be  more  robust  to  noise  than  the  one-dimensional  model 

In  the  experiments  of  the  previous  chapter  we  considered  each  densitometru 
profile  independently  of  the  adjacent  profiles,  i  e  the  spatial  continuitv  of  the  ves- 
"el  and  the  bac  kground  wu  not  taken  into  account  In  chapter  <  we  disc  us  sec  j  a 
heuristic  wav  to  use  the  spatial  conlinuitv  of  the  vessel  with  the  one  dimensional 
,.>p  roach  This  is  to  ht  splines  to  the  est  i  mated  v  esse  I  parameter  M-cpienc  e*  We  will 
refer  to  this  modified  approach  as  the  smoothed  one -dimensional  approa<  h  In  order 
to  <  ompare  it  t o  t  he  t  wo-dirnensional  approac  h  we  w  ill  use  t  he  same  spline  param 
eters  for  the  smoothed  one-dimensional  approac  h  as  th<e*e  of  the  t wo-dimensional 
approach  In  the  following  experiments  we  will  compare  the  t wo-dimensional  ap¬ 
proac  h  with  both  t he  one-dtmens»onal  approach  and  the  smoothed  one-dtmensionat 


approach. 

Similarly  to  the  previous  chapter,  the  algorithms  were  tested  on  computer  gen¬ 
erated  data  (synthetic  data),  on  x-rays  of  contrast-medium-filled  cylindrical  phan¬ 
toms.  and  on  real  coronary  angiograms. 


5.1  Synthetic  Data 

First  we  examine  the  performance  of  the  model  on  synthetic  data.  The  data  were 
generated  assuming  the  two-dimensional  model  is  valid.  We  conducted  a  number 
of  experiments  to  study  the  properties  of  the  two-dimensional  modeling  approach 
and  compare  it  with  the  one-dimensional  approach.  In  particular,  we  investigated 
the  behavior  of  the  two  approaches  for  different  levels  of  noise  and  different  vessel 
dimensions. 

In  order  to  be  able  to  make  quantitative  statements  about  the  experiments,  we 
t onsidered  first  vessels  of  constant  radius  (r  q  r  )  and  alpha  (a  q  -  o„).  Also, 
we  rfiose  the  axes  of  the  vessels  to  be  straight  vertical  lines  (C|  q  c„ .  c2  q  -  q) 
This  makes  it  possible  to  isolate  the  effects  that  the  level  of  the  noise,  the  radius 
of  the  vessel  and  the  parameter  alpha  have  on  the  problem  It  also  simplifies  the 
implerneritat  ion  of  t  fie  one-dimensional  algorithm,  and  makes  the  comparison  easier 
i  n<  e  the  two-dimensional  algorithm  does  not  take  these  facts  into  consideration, 
t  hev  should  not  bias  our  experiments  A  more  realistic  example  with  parameter 
sequences  obtained  from  a  real  angiogram  will  fie  considered  at  t  fie  end  of  this 
sec  t  ion 

lii  eac  fi  of  the  experiment'  we  desc  ribe  below  all  t  fie  parameter-  (background, 
c  r  ii  -  and  noise  variance;  were  fixed  except  one  tSe  considered  one  noise 

-ample  tor  e«n  h  value  of  the  variable  parameter  and  e-iimaied  the  vessel  parameter 
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Figure  -VI:  Two-Dimensional  Synthetic  Image 


and  hence  the  vessel  cross-sections  are  constant,  we  can  compute  the  average 
and  standard  deviation  of  the  radius  and  area  sequence  estimates  over  q  This 
is  similar  to  the  experiments  of  chapter  4,  where  we  considered  K()  noise  samples 
of  a  one-dimensional  profile  and  averaged  the  estimates  The  standard  deviation 
was  normalized,  that  is.  it  was  divided  by  the  true  value  of  the  parameter  being 
estimated  I  radius  or  area) 

V  typical  example  of  a  synthetic  test  image  is  shown  in  figure  VI  The  back¬ 
ground  is  a  fifth  degree  two-dimensional  polynomial,  obtained  from  a  real  an¬ 
giogram  The  vessel  parameter  values  are  r,  q  2~>  pixels,  r  q  h  pixels.  «  q  2. 
for  all  q.  and  >  diay{l.  I|  (pixels)*  The  standard  deviation  of  the  noise  is  a  2. 
w  h ii  h  i  orres ponds  to  signal  to  noise  ratio  of  1 f>  2  d M 

The  degree  of  bo!  h  t  fi«  two-dimensional  and  the  one  dimensional  polynomials 
used  ili  the  corresponding  estimation  algorithms  was  V  the  same  as  that  of  the 
svnthetic  cfata  I'fie  spac  mg  of  tfie  nodes  of  the  cubic  splines  used  in  the  last  stage 
of  the  two-dimensional  algorithm  ami  for  MiKMitliitiK  the-  one  dimensional  estimate's 


5.1.1  Effect  of  Noise  Level 


for  t  h i''  experiment  dll  the  parameters  were  fixed  except  the  noise  level  The  syn- 
thetn  image  we  used  was  tfiat  of  figure  5.1  ’’  he  vessel  parameters  were  ct  q  25 

pixels,  rq  *»  pixels,  and  o  q  2.  for  all  >,  and  >'  dt<i</(l.l)  (pixels)*  To 

this  profile  we  added  7  different  levels  of  noise.  For  each  noise  sample  we  estimated 
the  vessel  parameter  sequences  using  the  two-dimensional  and  one-dimensional  ap¬ 
proaches  The  signal  to  noise  ratio  (SNR)  is  defined  as  in  equation  (4.1)  of  the 
previous  chapter. 

The  average  of  the  radius  estimates  as  a  function  of  signal  to  noise  ratio  is  shown 
in  figure  5.2(a)  and  the  normalized  standard  deviation  is  shown  in  figure  5.2(b). 
The  solid  line  represents  the  two-dimensional  model  estin  s.  the  dot-dashed  line 
represents  the  one-dimensional  model  estimates  without  any  smoothing,  and  the 
dashed  line  represents  the  smoothed  one-dimensional  estimates.  The  true  radius 
is  shown  m  dotted  line.  We  observe  that  the  average  of  the  estimates  for  the 
t  Ao-dimensional  approach  is  very  similar  to  that  of  the  one-dimensional  approach 
I  f .«  smoothing  of  t  fie  one-dimensional  estimates  does  not  affect  the  average  (the 
dot  -dashed  and  dashed  lines  coincide),  but.  as  expected,  improves  considerably  the 
standard  deviation  I'he  standard  deviation  of  the  two-dimensional  estimates  is 
oniv  -light  I  v  better  than  that  of  t  he  smoothed  one-dimensional  estimates. 

I'he  average  of  t  fie  area  est  i mates  as  a  fum  t  ion  of  signal  to  noise  rat  10  is  shown 
,n  hgur»  ’>  !(a|  and  the  normalized  standard  deviation  is  shown  in  figure  >.d(b). 
I  hi  -olid  line  represents  t  fie  two-dimensional  model  estimates,  the  dot-dashed  line 
•  t  iow  -  the  orieofimensional  model  est  imat  os  w  it  limit  an  v  smoot  h  i  tig.  ami  t  fie  dashed 
Mu  ■■  f  ,ow  s  the  smoot  fied  one-dimensional  estimates  I  he  dotted  line  represents  the 
true  area  In  figure  '«  .Pal  the  dot-dashed  and  dashed  lines  toincide  I  h«*  two- 
dimensional  approa<  h  offers  some  improvement  in  t he  average  of  the  estimates, 
and  on  i  \  a  'tnall  improvement  in  t  fie  standard  deviation  of  t  he  smoothed  ofie- 
d  imerisiolial  est  itliate-- 
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Figure  5.2:  Radius  Estimation  for  Different  SNR's 
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Figure  5.3:  Area  Estimation  for  Different  SNR’s 
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5.1.2  Effect  of  Vessel  Radius 


For  this  experiment  all  the  parameters  were  fixed  except  the  vessel  radius.  The 
vessel  parameters  were  cx  q  =25  pixels,  r  q  =  r0  pixels,  and  a\q\  -  2,  for  all  q ,  and 
5  -  dtag(l,l)  (pixels)3.  The  standard  deviation  of  the  noise  was  fixed  at  <7  =  2  (this 
corresponds  to  signal  to  noise  ratio  of  16.2  dB  when  r0  -  8  pixels).  The  values  of 
the  radius  were  r0  =  3,5,8. 10  pixels.  We  estimated  the  vessel  parameter  sequences 
for  each  radius  using  the  two-dimensional  and  one-dimensional  approaches. 

The  average  of  the  radius  estimates  as  a  function  of  radius  is  shown  in  fig¬ 
ure  5.4(a)  and  the  normalized  standard  deviation  is  shown  in  figure  5.4(b).  The 
solid  line  represents  the  two-dimensional  model  estimates,  the  dot-dashed  line  repre¬ 
sents  the  one-dimensional  model  estimates  without  any  smoothing,  and  the  dashed 
line  represents  the  smoothed  one-dimensional  estimates.  The  true  radii  lie  on  the 
diagonal  (dotted  line).  We  observe  that  the  average  of  the  estimates  for  the  two- 
dimensional  approach  is  almost  identical  to  that  of  the  one-dimensional  approach, 
and  both  are  close  to  the  true  value.  The  two-dimensional  approach  offers  some 
improvement  in  the  standard  deviation  of  the  estimates,  especially  for  smaller  radii. 

The  average  of  the  area  estimates  as  a  function  of  area  is  shown  in  figure  5.5(a) 
and  the  normalized  standard  deviation  is  shown  in  figure  5.5(b).  Both  axes  of 
figure  5.5(a)  and  the  x-axis  of  figure  5.5(b)  show  the  square  root  of  the  area,  so 
that  the  data  points  are  more  evenly  spaced.  The  solid  line  represents  the  two- 
dimensional  model  estimates,  the  dot-dashed  line  shows  the  one-dimensional  model 
estimates  without  any  smoothing,  and  the  dashed  line  shows  the  smoothed  one¬ 
dimensional  estimates.  The  dotted  line  represents  the  true  area.  Again,  the  average 
of  the  estimates  for  the  two-dimensional  approach  is  very  close  to  that  of  the  one¬ 
dimensional  approach,  and  both  are  close  to  the  true  area.  The  two-dimensional 
approach  offers  some  improvement  in  the  standard  deviation  of  the  estimates  over 
the  smoothed  one-dimensional  approach. 
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Figure  5.4:  Radius  Estimation  for  Different  Radii 
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5.1.3  Effect  of  Parameter  Alpha 

For  this  experiment  all  the  parameters  were  fixed  except  a  The  vessel  parameters 
were  ct  q  25  pixels,  r  q  -  H  pixels,  and  a  q  o  .  for  all  q.  and  >  (1.1) 

(pixels*)  The  standard  deviation  of  the  noise  was  fixed  at  r  2  Ithis  corresponds 
to  signal  to  noise  ratio  of  16.2  dB  when  a  2)  The  parameter  alpha  ranged  from 
a  0.5.  to  q  2.5  in  steps  of  0.5.  For  each  a  we  estimated  the  vessel  parameter 
sequences  using  the  two-dimensional  and  one-dimensional  approaches 

The  average  of  the  radius  estimates  as  a  function  of  alpha  is  shown  in  fig¬ 
ure  5.6(a)  and  the  normalized  standard  deviation  is  shown  in  figure  5.6(b)  The 
solid  line  represents  the  two-dimensional  model  estimates,  the  dot-dashed  line  repre¬ 
sents  the  one-dimensional  model  estimates  without  any  smoothing,  and  the  dashed 
line  represents  the  smoothed  one-dimensional  estimates.  The  true  radius  is  shown 
in  dotted  line.  We  observe  that  the  average  of  the  estimates  for  the  two-dimensional 
approach  is  almost  identical  to  that  of  the  one-dimensional  approach,  and  both  are 
close  to  the  true  value.  The  two-dimensional  approach  offers  some  improvement  in 
the  standard  deviation  of  the  estimates  only  when  alpha  is  small. 

The  average  of  the  area  estimates  as  a  function  of  alpha  is  shown  in  figure  5.7(a) 
and  the  normalized  standard  deviation  is  shown  in  figure  5.7(b).  The  solid  line 
represents  the  two-dimensional  model  estimates,  the  dot-dashed  line  shows  the  one¬ 
dimensional  model  estimates  without  any  smoothing,  and  the  dashed  line  shows 
the  smoothed  one-dimensional  estimates.  The  dotted  line  represents  the  true  area. 
Again,  the  average  of  the  estimates  for  the  two-dimensional  approach  is  very  close 
to  that  of  the  one-dimensional  approach,  and  both  are  close  to  the  true  area.  The 
two-dimensional  approach  offers  some  improvement  in  the  standard  deviation  of  the 
estimates  over  the  smoothed  one-dimensional  estimates,  more  so  for  small  alpha. 
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Figure  5.6:  Radius  Estimation  for  Different  Alphas 
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Figure  5.8:  Two-Dimensional  Synthetic  Image 


5.1.4  Realistic  Example 

The  previous  experiments  were  designed  to  study  the  behavior  of  the  two  approaches 
for  different  levels  of  noise  and  different  vessel  dimensions.  Here  we  present  a  more 
realistic  example  in  order  to  compare  the  overall  performance  of  the  two  methods 
on  synthetic  data. 

The  synthetic  image  we  used  in  this  experiment  is  100  x  100  pixels,  and  is  shown 
in  figure  5.8.  No  noise  has  been  added  to  this  image.  The  background  is  a  fifth 
degree  two-dimensional  polynomial,  obtained  from  a  real  angiogram  profile.  The 
vessel  parameter  sequences  Cilg),  r[g],  and  ajg]  are  shown  in  figures  5.9(a),  (b),  and 
(c),  respectively.  They  too  were  chosen  to  closely  match  the  vessel  projection  of  a 
real  angiogram.  The  covariance  matrix  of  the  blurring  function  is  S  =  diag(l,l) 
(pixels)2.  The  units  of  q  and  the  sequences  ci [qrj  and  r[g]  are  in  pixels. 

To  the  synthetic  image  of  figure  5.8  we  added  seven  different  levels  of  noise.  The 
parameter  sequences  of  a  section  of  this  image  (80  x  80  pixels,  shown  in  dotted  lines 


in  figure  5.9)  were  estimated  using  the  two  approaches.  The  degree  of  both  the 
two-dimensional  and  the  one-dimensional  polynomials  used  in  the  corresponding 
approaches  was  5,  the  same  as  that  of  the  synthetic  data.  The  spacing  of  the  nodes 
of  the  cubic  splines  used  in  the  last  stage  of  the  two-dimensional  algorithm  and  for 
smoothing  the  one-dimensional  estimates  was  7  points.  A  second  iteration  of  the 
one-dimensional  approach  was  performed  to  get  perpendicular  profile  estimates.  A 
typical  example  of  the  vessel  parameter  estimates  is  shown  in  figure  5.10.  They 
correspond  to  the  case  of  signal  to  noise  ratio  of  16.3  dB.  The  vessel  parameter 
sequences  ci[g],  r[q],  and  aj^]  are  shown  in  figures  5.10(a),  (b),  and  (c),  respectively. 
The  solid  line  represents  the  two-dimensional  model  estimates,  the  dot-dashed  line 
represents  the  one-dimensional  model  estimates,  and  the  dotted  line  shows  the  true 
parameter  sequences.  The  smoothed  one-dimensional  estimates  are  not  shown  in 
this  figure. 

The  differences  between  the  true  and  estimated  radius  and  area  sequences  were 
obtained.  We  then  computed  the  average  and  standard  deviation  of  these  differ¬ 
ence  sequences.  The  standard  deviation  was  normalized  by  the  average  value  of 
the  true  radius  and  area  sequences,  respectively.  The  results  are  shown  in  fig¬ 
ures  5.11  and  5.12  as  a  function  of  signal  to  noise  ratio.  The  signal  to  noise 
ratio  (SNR)  is  defined  as  in  equation  (4.1)  of  the  previous  chapter.  The  solid  line 
represents  the  two-dimensional  model  estimates,  the  dot-dashed  line  represents  the 
one-dimensional  model  estimates  without  any  smoothing,  and  the  dashed  line  repre¬ 
sents  the  smoothed  one-dimensional  estimates.  Since  we  are  plotting  the  difference 
between  the  true  and  estimated  values,  the  perfect  result  would  be  the  zero  se¬ 
quence  shown  in  dotted  line.  For  both  the  radius  and  the  area  estimates  we  observe 
that  the  two-dimensional  approach  offers  some  improvement  in  the  average  and  the 
standard  deviation  of  the  estimates. 

The  results  of  this  more  realistic  example  are  similar  to  those  of  the  previous 
sections.  The  difference  in  performance  between  the  two-dimensional  approach  and 


Figure  5.10:  Vessel  Parameter  Estimates 
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the  smoothed  one-dimensional  approach  we  observe  here  is  comparable  to  that  of 
the  previous  sections.  This  means  that  the  earlier  experiments  were  reasonable, 
even  though  the  vessel  shape  was  too  simple. 


5.1.5  Remarks 


Summarizing  the  results  of  the  synthetic  data  experiments,  we  note  that  the  two- 
dimensional  approach  offers  some  improvement  in  the  variance  of  the  estimates, 
while  the  average  is  very  close  to  the  true  value  for  both  approaches.  As  expected, 
the  smoothing  of  the  one-dimensional  estimates  reduces  the  standard  deviation 
considerably,  but  does  not  affect  the  average. 


5.2  Vessel  Phantoms 


Now  we  examine  the  performance  of  the  algorithms  on  the  vessel  phantom  x-rays. 
The  phantoms  are  the  same  ones  that  we  used  in  chapter  4.  For  each  vessel  phantom 
we  considered  a  section  of  the  x-ray  which  contained  50  consecutive  profiles  of  the 
phantom,  and  estimated  the  vessel  parameter  sequences  using  the  one-dimensional 
and  two-dimensional  approaches.  Since  the  phantoms  have  constant  cross-sections, 
we  can  compute  the  average  and  the  standard  deviation  of  the  radius  and  area 
estimates  over  q.  The  standard  deviation  of  the  estimates  was  normalized  in  the 
same  way  as  in  section  4.2,  by  the  values  on  a  diagonal  line  close  to  the  average 
of  the  one-dimensional  estimates  (the  exact  correspondence  between  the  radius  in 
millimeters  and  the  radius  in  pixels  is  not  known).  This  experiment  is  similar  to 
the  radius  experiment  in  the  synthetic  data  case  (section  5.1.2).  The  degree  of  the 
polynomials  used  in  the  estimation  algorithms  was  3.  As  we  discussed  in  section  4.2, 
a  third  degree  polynomial  to  model  an  approximately  flat  background  has  better 
performance.  The  spacing  of  the  nodes  of  the  cubic  splines  used  in  the  last  stage 
of  the  two-dimensional  algorithm  and  for  smoothing  the  one-dimensional  estimates 
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was  7  points. 

The  average  of  the  radius  estimates  as  a  function  of  true  radius  is  shown  in 
figure  5.13(a)  and  the  normalized  standard  deviation  is  shown  in  figure  5.13(b). 
The  solid  line  represents  the  two-dimensional  model  estimates,  the  dot-dashed  line 
represents  the  one-dimensional  modei  estimates  without  any  smoothing,  and  the 
dashed  line  represents  the  smoothed  one-dimensional  estimates.  The  units  of  the 
x-axis  are  millimeters  (mm)  and  of  the  y-axis  pixels.  As  we  discussed  in  the  pre¬ 
vious  chapter,  we  are  interested  in  relative  measurements.  Thus  the  estimates  are 
satisfactory  when  they  lie  on  a  straight  line  passing  through  the  origin.  We  observe 
that  the  average  of  the  radius  estimates  for  the  two-dimensional  model  roughly  lies 
on  a  straight  line  through  the  origin,  and  is  very  close  to  the  one-dimensional  model 
average.  The  standard  deviation  for  the  two-dimensional  approach  is  lower  than 
the  standard  deviation  for  the  smoothed  one-dimensional  approach. 

The  average  of  the  area  estimates  as  a  function  of  area  is  shown  in  figure  5.14(a) 
and  the  standard  deviation  is  shown  in  figure  5.14(b).  Both  axes  of  figure  5.14(a) 
and  the  x-axis  of  figure  5.14(b)  show  the  square  root  of  the  area,  so  that  the  data 
points  are  more  evenly  spaced.  The  solid  line  represents  the  two-dimensional  model 
estimates,  the  dot-dashed  line  shows  the  one-dimensional  model  estimates  without 
any  smoothing,  and  the  dashed  line  shows  the  smoothed  one-dimensional  estimates. 
Again,  there  is  very  little  difference  in  the  average  of  the  estimates  of  the  two 
approaches.  The  behavior  of  the  average  area  estimate  was  discussed  in  section  4.2 
of  the  previous  chapter.  The  standard  deviation  of  the  two-dimensional  model 
estimates  is  lower  than  the  standard  deviation  of  the  smoothed  one-dimensional 
model  over  most  of  the  range  of  area  values. 
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(b)  Normalized  Standard  Deviation  of  Radius  Estimates 


Figure  5.13:  Phantom  Experiment:  Radius  Estimation 
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Figure  5.14:  Phantom  Experiment:  Area  Estimation 


5.3  Coronary  Angiograms 


The  performance  of  the  algorithms  was  also  tested  on  real  coronary  angiograms. 
As  we  mentioned  earlier,  it  is  very  difficult  to  evaluate  the  performance  of  the 
algorithms  in  this  case,  because  we  do  not  know  the  true  dimensions  of  the  arteries. 
However,  a  simple  inspection  of  the  angiograms  and  their  models  can  give  us  a  good 
indication  of  the  relative  performance  of  the  two  approaches.  Thus  the  comparisons 
of  this  section  will  be  qualitative,  in  contrast  to  the  quantitative  results  of  the 
previous  sections. 

In  the  examples  we  present,  the  approximate  direction  of  the  blood  vessels  is 
vertical.  The  one-dimensional  algorithm  was  applied  to  the  profiles  perpendicular 
to  this  direction.  More  iterations  of  the  algorithm  provided  the  parameters  of 
the  perpendicular  vessel  cross-sections.  The  degree  of  the  polynomials  used  in 
the  estimation  algorithms  was  5.  The  spacing  of  the  nodes  of  the  cubic  splines 
used  in  the  last  stage  of  the  two-dimensional  algorithm  and  for  smoothing  the  one¬ 
dimensional  estimates  was  7  points. 

The  performance  of  the  two-dimensional  model  is  shown  in  figure  5.15.  Part  (a) 
shows  the  original  image,  which  is  80  x  80  pixels;  part  (b)  shows  the  model;  part 
(c)  shows  the  vessel  component  of  the  model;  and  part  (d)  shows  the  background 
component.  The  estimated  radius  and  area  sequences  are  shown  in  solid  line  in 
figures  5.16(a)  and  5.16(b)  respectively. 

The  performance  of  the  one-dimensional  model  is  shown  in  figure  5.17.  Part 
(a)  shows  the  original  image;  part  (b)  shows  the  model;  parts  (c)  and  (d)  show  the 
vessel  and  background  components  of  the  model,  respectively.  Figure  5.18  shows 
the  smoothed  one-dimensional  model.  In  these  figures  we  show  the  results  of  the 
first  iteration  of  the  algorithm.  The  profiles  in  the  first  stage  of  the  optimization 
are  parallel  and  nonoverlapping,  and  hence  can  be  displayed  in  an  image.  The 
profiles  of  the  second  stage  are  in  different  directions  and  can  only  be  plotted  one 
at  a  time.  The  displayed  results  are  adequate  for  a  qualitative  comparison  of  the 
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Figure  5.15:  Two-Dimensional  Model  (Kxample  A) 
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Figure  5.16:  Estimated  Parameter  Sequences  (Example  A) 
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Figure  5.18:  Smoothed  One-Dimensional  Model  (Example  A) 
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algorithm  performances.  For  the  one-dimensional  model,  the  estimated  radius  and 
area  sequences  of  the  profiles  perpendicular  to  the  vessel  are  shown  in  dot-dashed 
line  in  figures  5.16(a)  and  5.16(b),  respectively.  The  smoothed  radius  and  area 
sequences  are  shown  in  dashed  lines. 

Comparing  figures  5.15  and  5.18,  we  observe  that  both  models  look  pretty  close 
to  the  original.  As  should  be  expected  the  background  in  the  one-dimensional 
model  is  discontinuous  in  the  direction  of  the  vessel.  However,  when  we  look  at  the 
vessel  component  of  the  one-dimensional  model,  we  see  a  valley  in  the  intensity  at 
the  narrow  point  in  the  middle.  This  valley  is  not  present  in  the  two-dimensional 
vessel  model,  and  is  apparently  caused  by  the  small  vessels  in  the  background.  This 
indicates  that  the  two-dimensional  model  is  more  robust  to  non-white  noise  in  the 
background. 

Another  example  is  shown  in  figures  5.19,  5.20,  and  5.21;  the  resolution  is  80x80 
pixels.  We  observe  that  the  vessel  component  of  the  one-dimensional  model  has 
some  wide  points  which  are  obviously  caused  by  the  background  noise.  The  effect  of 
the  noise  on  the  two-dimensional  vessel  estimate  is  considerably  smaller.  The  signal 
to  noise  ratio  is  very  low  in  this  case,  and  the  performance  of  the  two-dimensional 
model  is  clearly  superior.  This  observation  is  consistent  with  the  results  of  the 
previous  sections,  where  we  showed  that  the  standard  deviation  of  the  estimates 
in  the  two-dimensional  approach  becomes  significantly  smaller  for  small  signal  to 
noise  ratios. 

The  estimated  radius  and  area  sequences  of  the  profiles  perpendicular  to  the 
vessel  are  shown  in  figures  5.22(a)  and  5.22(b),  respectively.  The  two-dimensional 
model  estimates  are  in  solid  lines,  the  one-dimensional  model  estimates  are  in 
dot-dashed  lines,  and  the  smoothed  one-dimensional  model  estimates  are  shown 
in  dashed  lines. 
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Figure  5.20:  One-Dimensional  Model  (Example  B) 
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Figure  5.21:  Smoothed  One- Dimensional  Model  (Example  B) 
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5.4  Summary 


The  synthetic  and  phantom  data  experiments  we  presented  in  this  chapter  in¬ 
dicate  that  when  compared  to  the  smoothed  one-dimensional  approach  the  two- 
dimensional  approach  offers  some  improvement  in  the  variance  of  the  estimates, 
while  the  average  is  very  close  to  the  true  value  for  both  approaches.  When  we 
examined  the  behavior  of  the  algorithms  on  real  coronary  angiograms,  however, 
we  found  that  the  two-dimensional  model  is  more  robust  to  non-white  noise  in  the 
background. 

As  we  discussed  in  chapter  3,  the  two-dimensional  approach  requires  a  lot  more 
computation  than  the  one-dimensional  approach.  Again,  when  one  considers  the 
time  and  effort  that  is  required  to  get  an  angiogram  of  the  coronary  arteries,  it 
appears  to  be  reasonable  to  spend  more  time  to  get  better  estimates  of  the  severity 
of  the  obstructions. 


Chapter  6 


Summary  and  Future  Work 

6.1  Summary 

In  this  thesis  we  presented  a  modeling  approach  for  the  estimation  of  the  dimen¬ 
sions  of  the  coronary  arteries  from  coronary  angiograms.  The  model  we  developed 
accounts  for  the  structure  of  the  background,  as  well  as  distortions  introduced  by 
the  imaging  system.  The  model  takes  into  account  the  continuity  of  both  the  vessel 
and  the  background.  The  desired  vessel  dimensions,  that  is,  the  diameter  and  cross- 
sectional  area  at  each  point  along  the  vessel,  can  be  obtained  easily  from  the  model 
parameters.  Two  models  were  presented.  The  first  is  a  two-dimensional  model  of  a 
section  of  the  angiogram  containing  a  single  vessel.  The  second  is  one-dimensional 
model  of  the  profiles  of  the  lines  crossing  the  vessel. 

The  maximum  likelihood  estimates  of  the  model  parameters  were  obtained  by 
iterative  algorithms.  The  computational  requirements  for  the  one-dimensional  al¬ 
gorithm  are  reasonable.  The  two-dimensional  algorithm  requires  much  more  com¬ 
putation  than  the  one-dimensional  algorithm.  The  amount  of  computation  for  both 
algorithms  is  considerably  greater  than  that  of  the  slope  methods. 

We  tested  the  algorithms  on  three  types  of  data.  The  first  type  is  computer  gen¬ 
erated  data.  The  second  is  x-rays  of  contrast-medium-fiiled  cylindrical  phantoms. 


The  algorithms  were  also  tested  on  real  coronary  angiograms.  In  the  experiments 
we  presented  we  first  compared  the  performance  of  the  one-dimensional  model  to 
the  slope  approaches.  Then  we  compared  the  performance  of  the  one-dimensional 
model  to  that  of  the  two-dimensional  model. 

Our  results  indicate  that  both  the  one-dimensional  and  the  two-dimensional 
algorithm  have  better  performance  than  the  conventional  slope  approaches.  This 
is  because  the  modeling  approach  uses  ail  the  information  in  the  projection  image 
to  estimate  the  model  parameters,  while  the  slope  approaches  use  only  the  slope 
at  the  sides  of  the  vessel  to  determine  the  vessel  boundaries  and  to  estimate  the 
background.  We  also  showed  that  the  two-dimensional  algorithm  is  more  robust  to 
non-white  noise  in  the  background  than  the  one-dimensional  algorithm. 

As  we  discussed  in  the  thesis,  the  cross-sections  of  the  coronary  obstructions 
are  frequently  irregular  in  shape.  This  suggests  that  more  complicated  models  of 
the  vessel  would  be  desirable.  As  the  model  gets  more  complicated,  however,  the 
problem  of  estimating  the  parameters  from  one  projection  becomes  more  and  more 
ill-conditioned.  Also  the  amount  of  computation  increases.  The  model  we  developed 
in  this  thesis  is  an  attempt  to  balance  the  need  for  an  accurate  model  of  the  vessels 
and  the  requirement  of  a  well-conditioned  estimation  problem. 

Our  model  can  easily  be  extended  to  handle  the  case  of  multiple  vessels  in  the 
image  segment.  It  can  also  be  extended  to  use  more  projections  to  reconstruct 
the  elliptical  cross-sections.  The  problem  in  both  cases  of  course  is  that  the  com¬ 
putational  requirements  and  the  complexity  of  the  estimation  procedure  increase 
significantly. 

6.2  Suggestions  for  Further  Research 

A  problem  that  is  not  treated  in  this  thesis  is  that  of  estimating  the  severity  of 
obstructions  that  occur  at  points  where  the  vessel  separates  into  two  branches. 


119 


This  is  an  important  problem,  because  stenoses  can  occur  at  vessel  bifurcations  |4  . 
To  estimate  the  artery  dimensions  at  a  bifurcation,  one  could  extend  the  modeling 
approach  of  this  thesis.  The  vessel  component  of  the  model  would  have  to  be 
modified  to  describe  a  bifurcation.  However,  the  new  model  will  most  likely  involve 
more  parameters,  making  the  parameter  estimation  procedure  more  complicated. 

Another  interesting  problem  is  that  of  reconstructing  the  vessel  cross-sections 
from  two  or  more  projections.  We  have  already  discussed  the  problems  with  ob¬ 
taining  more  than  one  view  of  the  coronary  arteries  and  the  difficulty  of  proper 
registration  of  the  projections  in  the  various  views.  Despite  these  problems  the 
problem  of  vessel  reconstruction  from  two  or  more  projections  has  been  of  interest 
to  many  researchers.  According  to  the  projection-slice  theorem,  if  we  know  the 
projections  at  all  angles,  we  can  reconstruct  the  object  exactly.  However,  the  vessel 
reconstruction  problem  has  a  very  special  structure.  Thus  it  might  be  possible  to 
reconstruct  the  cross-section  from  two  orthogonal  projections.  The  conditions  under 
which  this  can  be  achieved  are  not  known.  The  general  problem  of  reconstructing 
an  array  of  zeros  and  ones  from  two  projections  (the  sum  of  its  rows  and  columns) 
is  considered  in  |7j.  Approximate  reconstructions  have  been  considered  by  Reiber 
et  a  I.  18.33  and  Spears  et  a  I.  [33] .  Spears  et  a/,  use  a  maximum  entropy  iterative 
algorithm  to  reconstruct  the  cross-sections  from  a  small  number  of  views.  All  of 
the  above  approaches  assume  that  the  vessel  component  of  each  projection  can  be 
isolated  from  the  background  in  the  presence  of  blurring  and  noise.  However,  as  we 
argued  in  this  thesis,  the  problem  of  isolating  the  vessel  component  of  the  projection 
is  not  easy. 

Finally,  there  is  the  problem  of  combining  the  information  in  consecutive  frames 
in  time.  This  would  require  an  accurate  procedure  to  estimate  the  movement  of  the 
vessel  from  frame  to  frame,  and  could  lead  to  improved  parameter  estimates. 


Bibliography 


lj  H.  L.  Abrams,  ed.  Coronary  Arteriography.  Little,  Brown  and  Company, 
1983. 

2'  G.  J.  Agin,  T.  O.  Binford,  “Computer  Description  of  Curved  Objects,”  IEEE 
Trans.  Comp.,  Vol.  C-25,  No.  4,  pp.  439-449,  April  1976. 

3!  E.  L.  Alderman  et  al.,  “Quantitation  of  Coronary  Artery  Dimensions  Using 
Digital  Image  Processing,”  Digital  Radiography,  SPIE  Vol.  314,  pp.  273-278, 
1981. 

4;  M.  Amiel  et  al.,  Coronary  Artery  Diseases.  Springer- Verlag,  1984. 

5;  K.  Barth  et  al.,  “The  Objective  Measurement  of  Coronary  Obstructions  by 
Digital  Image  Processing,”  Proc.  of  Int.  Workshop  on  Physics  and  Eng.  in 
Med.  Imaging,  IEEE  Cat.  No.  82  CH1751-7,  pp.  74-78,  1982. 

6;  K.  Barth  et  al.,  “Quantification  of  Coronary  Stenoses,”  in  Radiological  Func¬ 
tional  Analysis  of  the  Vascular  System ,  (Ed.  F.  H.  W.  Heuck),  Springer- Verlag, 
Berlin,  pp.  207-213,  1983. 

7'  S.  K.  Chang,  Y.  R.  Wang,  “Three-Dimensional  Object  Reconstruction  from 
Orthogonal  Projections,”  Pattern  Recognition,  Vol.  7,  pp.  167-176,  1975. 

8  S.  D.  Conte,  Carl  De  Boor,  Elementary  Numerical  Analysis.  McGraw-Hill, 
1980. 

9  D.  W.  Crawford  et  al.,  “Rotating  Step-Wedge  Technique  for  Extraction  of 
Luminal  Cross-Sectional  Area  Information  from  Single  Plane  Coronary  An¬ 
giograms,”  Investigative  Radiology,  Vol.  12,  pp.  307-313,  1977. 

10;  Carl  De  Boor,  A  Practical  Guide  to  Splines.  Springer- Verlag,  Applied  Math¬ 
ematical  Sciences,  Vol  27,  1978. 

llj  T.  A.  De  Rouen  et  al.,  “Variability  in  the  Analysis  of  Coronary  Arteriograms,” 
Circulation,  55,  pp.  324-328,  1977. 


121 


12  K.  M.  Detre  et  a I.,  “Observer  Agreement  in  Evaluating  Coronary  Angiograms,” 
Circulation,  52,  pp.  979-986.  1975. 

13  K.  L.  Gould,  K.  Lipscomb,  “Effects  of  Coronary  Stenoses  on  Coronary  Flow 
Reserve  and  Resistance,"  Am.  J.  of  Cardiology,  Vol.  34,  pp.  48-55,  July  1974. 

14  K.  L.  Gould.  "Pressure- Flow  Characteristics  of  Coronary  Stenoses  in  Unse¬ 
dated  Dogs  at  Rest  and  During  Coronary  Vasodilation,"  Circ.  Res.,  VoJ.  43, 
pp.  242-..,  1978. 

15  C.  J.  Kooijman  et  a  I.,  "Computer-Aided  Quantitation  of  the  Severity  of  Coro¬ 
nary  Obstructions  from  Single  View  Cineangiograms,"  Proceedings  of  Int. 
Symp.  on  Medical  Imaging  and  Image  Interpretation,  pp.  59-64,  1982. 

16  J.  H.  C.  Reiber  et  a/..  "Objective  Characterization  of  Coronary  Obstructions 
from  Cineangiograms  with  Single  and  Multiple  Views,"  Proceedings  of  Second 
Int.  Conf.  on  Visual  Psychophysics  and  Medical  Imaging,  Brussels,  Belgium, 
pp.  11-17,  1981. 

17  J.  H.  C.  Reiber  et  a/.,  "Transfer  Functions  of  the  X-Ray-Cine-Video  Chain 
Applied  to  Digital  Processing  of  Coronary  Cineangiograms."  in  Digital  Imaging 
in  Cardiovascular  Radiology,  (Eds.  P.  H.  Heintzen  and  R.  Brennecke).  Georg 
Thieme  Verlag,  Stuttgart,  pp.  89-104.  1983. 

18:  J.  F  C.  Reiber  et  a/.,  “3-D  Reconstruction  of  Coronary  Arterial  Segments  from 
Two  Projections,”  in  Digital  Imaging  in  Cardiovascular  Radiology,  (Eds.  P.  H. 
Heintzen  and  R.  Brennecke),  Georg  Thieme  Verlag,  Stuttgart,  pp.  151  163. 
1983. 

19  J.  H.  C.  Reiber  et  a/.,  “Coronary  Artery  Dimensions  from  Cineangiograms 
-  Methodology  and  Validation  of  a  Computer- Assisted  Analysis  Procedure." 
IEEE  Trans.  Med.  fmag.,  Vol.  MI-3,  n.  3,  pp.  131-140.  1984. 

20  J.  H.  C.  Reiber  et  al.,  “Computer  Assisted  Analysis  of  the  Severity  of  Obstruc¬ 
tions  from  Coronary  Cineangiograms:  A  Methodological  Review."  Automed- 
ic a,  Vol.  5,  pp.  219-238,  1984. 

21  J.  R.  Rice,  The  Approximation  of  Functions,  Vol  2.  Addison- Wesley,  1969. 

22  K.  Rossmann,  “Image-Forming  Quality  of  Radiographic  Screen-Film  Systems: 
The  Line  Spread-Function,”  Am.  J.  Roentgenology,  Rad.  Therapy  L  Nuclear 
Med..  Vol.  90,  No.  1,  pp.  886-906,  July  1963. 

23;  T.  Sandor  et  al.,  “Cine-Densitometric  Measurement  of  Coronary  Arterial 
Stenoses,"  Catheterization  and  Cardiovascular  Diagnosis,  5,  pp.  229-245,  1979. 


[24]  M.  E.  Sanmarco  et  a “Reproducibility  of  a  Consensus  Panel  in  the  Interpre¬ 
tation  of  Coronary  Angiograms,”  American  Heart  Journal,  96,  pp.  430-437, 
1978. 

[25]  H.  J.  Scudder,  “Introduction  to  Computer  Aided  Tomography,”  Proceedings 
IEEE,  Vol.  66,  n.  6,  pp.  628-637,  1978. 

[26]  R.  H.  Selzer  et  a ].,  “Computer  Analysis  of  Cardiovascular  Imagery,”  Proc.  of 
the  CALTEC/JPL  Conference  on  Image  Processing  Technology,  Data  Sources 
and  Software  for  Commercial  and  Scientific  Applications,  Pasadena,  California, 
pp.  6-1-6-20,  1976. 

[27]  K.  Shmueli,  “Estimation  of  Blood  Vessel  Boundaries  in  X-Ray  Images,”  Ph.D. 
Thesis,  Stanford  University,  1981. 

[28]  K.  Shmueli  et  a/.,  “Estimation  of  Blood  Vessel  Boundaries  in  X-Ray  Images,” 
Optical  Engineering,  Vol.  22,  No.  1,  pp.  110-116,  1983. 

29]  D.  N.  Smith  et  a  1.,  “A  Semiautomatic  Computer  Technique  for  Processing 
Coronary  Angiograms,”  Computers  in  Cardiology,  pp.  325-328,  1981. 

30!  J.  R.  Spears,  “Rotating  Step- Wedge  Technique  for  Extraction  of  Luminal 
Cross-Sectional  Area  Information  from  Single  Plane  Coronary  Angiograms,” 
Acta  Radiologica  Diagnosis ,  22,  pp.  217-225,  1981. 

31'  J.  R.  Spears  et  a I.,  “Computer  Aided  Densitometric  Evaluation  of  Coronary 
Cineangiograms,”  in  Radiological  Functional  Analysis  of  the  Vascular  System, 
(Ed.  F.  H.  W.  Heuck),  Springer- Verlag,  Berlin,  pp.  195-206,  1983. 

32:  J.  R.  Spears  et  ah,  “Computerized  Image  Analysis  for  Quantitative  Measure¬ 
ment  of  Vessel  Diameter  from  Cineangiograms,”  Circulation,  68,  pp.  453-461, 
1983. 

[33]  J.  R.  Spears  et  ah,  “Computer  Reconstruction  of  Luminal  Cross-Sectional 
Shape  from  Multiple  Cineangiographic  Views,”  IEEE  Trans.  Medical  Imag., 
Vol.  MI-2,  No.  1,  pp.  49-54,  March  1983. 

[34]  Harry  L.  Van  Trees,  Detection,  Estimation,  and  Modulation  Theory,  Part  I. 
Wiley,  New  York,  1968. 

[35]  E.  H.  Timothy  Whitten,  “Orthogonal  Polynomial  Trend  Surfaces  for  Irregularly 
Spaced  Data,”  Mathematical  Geology,  Vol.  2,  No.  2,  pp.  141-152,  1970. 

36]  L.  M.  Zir  et  ah,  “Interobserver  Variability  in  Coronary  Angiography,”  Circu¬ 
lation,  53,  pp.  627-632,  1976. 


123 


DISTRIBUTION  LIST 


i 

N 

§ 

W  Director 

aa  Defense  Advanced  Research  Project  Agency 

1400  Wilson  Boulevard 
•,*>  Arlington,  Virginia  22209 

jQ  Attn:  Program  Management 


£ 


Head  Mathematical  Sciences  Division 
Office  of  Naval  Research 
800  North  Quincy  Street 
Arlington,  Virginia  22217 


£ 


w 


Administrative  Contracting  Officer 
E19-628 

Massachusetts  Institute  of  Technology 
Cambridge,  Massachusetts  02  139 


Director 

Naval  Research  Laboratory 
Attn:  Code  2627 
Washington,  D.  C.  20375 


!! 

$ 


& 


Defense  Technical  Information  Center 
Bldg  5,  Cameron  Station 
Alexandria,  Virginia  22314 


Dr.  Judith  Daly 
DARPA  /  TTO 
1400  Wilson  Boulevard 
Arlington,  Virginia  22209 


Jfi 


£ 


DO  DA  AD  Code 


HX1241  (1) 


N00014  (1) 


N66017  (1) 


N00173  (6) 


S47031  (12) 


(1) 


S 


;  I*,  i; 


